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Abstract 


This  thesis  develops  a  method  to  analyze  the  maneuvering  forces  on  surfaced  and 
underwater  vehicles  with  complex  propulsors.  The  analysis  method  is  developed  for 
general  propellers  yet  has  unique  applicability  to  model  highly  contracting  stern  flows 
associated  with  integrated  propulsors.  Integrated  propulsors  exhibit  strong  coupling 
of  the  various  blade-rows  and  duct,  if  present,  to  the  vehicle  stern.  The  method 
developed  herein  provides  a  robust  means  to  analyze  propulsor-induced  maneuvering 
forces  including  those  arising  from  wake-adapted,  multi-stage,  ducted  propulsors. 

The  heart  of  the  maneuvering  force  prediction  is  a  three-dimensional,  unsteady 
lifting-surface  method  developed  as  the  first  part  of  this  thesis.  The  new  method  is 
designated  PUF-14  for  Propeller  Unsteady  Forces.  The  lifting-surface  method  uses 
many  advanced  techniques.  One  significant  advance  is  the  use  of  a  wake-adapted 
lattice  to  model  the  flow  through  the  propulsor.  In  related  research,  a  2-D  Kutta 
condition  has  been  augmented  using  Lagrangian  interpolation  to  dramatically  reduce 
the  required  computational  time  to  model  a  2-D  gust. 

The  second  thrust  of  this  thesis  couples  the  unsteady  lifting-surface  method  with  a 
three-dimensional,  time-average  Reynolds-Averaged  Navier-Stokes  flow  solver.  Rotat- 
ing a  propeller  through  a  spatially-varying  flow  field  causes  temporally-varying  forces 
on  the  propeller.  From  the  converged-coupled  solution,  the  maneuvering  and  blade- 
rate  forces  can  be  estimated.  This  thesis  explores  the  relationship  of  time-varying 
and  time-average  forces  in  the  flow  solver  and  potential-flow  domains.  Similarly,  it 
explores  the  relationship  of  the  effective  inflow  in  the  two  domains.  Finally,  this  the- 
sis details  the  synergistic  means  to  correctly  couple  the  potential-flow  method  to  a 
viscous  solver. 

Verification  and  validation  of  the  method  have  been  done  on  a  variety  of  geometries 
and  vehicles.  Preliminary  results  show  good  correlation  with  experiment.  The  results 
strongly  suggest  this  maneuvering  force  prediction  method  has  great  potential  for  the 
modern  propulsor  designer. 

Thesis  Supervisor:  Justin  E.  Kerwin 
Title:  Professor  of  Naval  Architecture 
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Chapter  1 

Introduction 

1.1  Overview 

In  order  to  study  complex  propulsors,  a  method  of  solving  the  resulting  time-dependent, 
non-linear  boundary  value  problem  is  required.   Methods  of  predicting  the  propeller 
blade-row  performance  in  steady  and  unsteady  flow  have  been  in  existence  for  many 
years.    Comprehensive  reviews  of  the  steady  and  unsteady  lifting-surface  theory  are 
given  in  Kerwin  [27]  and  Schwanecke  [34],  respectively. 

The  present  work  builds  on  advances  in  steady  lifting-surface  theory,  as  in  PBD- 
14  [28],  and  advances  in  unsteady  lifting-surface  theory,  such  as  PUF-2  [31].  By 
blending  the  best  features  of  the  two  methodologies,  a  robust  methodology  is  devel- 
oped to  calculate  the  time- varying  forces  on  highly  complex  propulsor  geometries. 
Then,  the  new  lifting-surface  methodology  is  mated  with  a  three-dimensional,  steady 
Reynolds- Averaged  Navier-Stokes  (RANS)  solver.  The  RANS  solver  captures  the 
viscous  nature  of  the  propulsor  characteristics  and  its  interaction  on  the  body.  The 
resulting  coupled  methodology,  documented  in  [42],  is  able  to  analyze  the  steady  and 
unsteady  blade-row  forces  on  highly  complex  propulsors. 

1.2  Need  for  a  New  Methodology 

Figures  1-1  and  1-2  illustrate  examples  of  highly  complex  propulsors.  The  hydrody- 
namic  complexity  arises  due  to  strong  interaction  of  the  rotor  with  the  full  afterbody. 
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The  duct,  if  present,  tends  to  inhibit  flow  separation  as  the  body  contracts  at  the 
vehicle  stern.  The  propulsor  components  are  strongly  coupled  and  must  be  designed 
and  analyzed  while  taking  into  account  their  full  interaction.  This  type  of  propulsor 
-  where  strong  interaction  occurs  between  components  and  the  body  -  is  known  as  an 
integrated  propulsor.  The  behavior  of  an  integrated  propulsor  during  a  maneuver  is 
difficult  to  predict  with  cylindrical-streamtube  design  methods.  The  current  research 
should  provide  a  relatively  fast  and  robust  analysis  method  for  integrated  propulsors. 
Figure  1-3  shows  the  typical  body  dynamics  during  a  maneuver.  During  a  given 
phase  of  the  turn,  the  body  sees  an  inflow  that  varies  slowly  compared  to  the  propeller 
rotation.  This  inflow  can  be  represented  with  a  time-average  RANS  solution  and 
the  propulsor  interaction  modeled  using  the  new  methodology.  Figures  1-1  and  1- 
5  illustrate  the  body-centered  coordinate  system  used  when  describing  maneuvering 
forces  and  moments. 


Figure  1-1:  Notional  submerged-body  integrated  propulsor  (Developed  at  MIT). 


Methodologies  to  calculate  unsteady  forces  on  highly  complex  propulsors,  such  as 
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Figure  1-2:  Notional  surface  ship  underwater  hull  form  with  an  integrated  propulsor  concept  (De- 
veloped at  NSWC  under  the  electric  drive  program.) 


figure  1-1,  exist  today  only  as  rudimentary  approximations.  The  new  methodology 
provides  a  method  to  assess  the  unsteady  forces  arising  in  the  wake-deficient  regions 
behind  control  surfaces,  etc.  and  those  arising  due  to  once-per-revolution  spatial 
variations  as  the  vehicle  maneuvers.  In  the  end,  the  coupled  methodology  should  be 
invaluable  to  the  modern  propulsor  designer. 

1.3     Functional  Description  of  the  Thesis 

The  areas  of  concentration  of  this  thesis  work  are  divided  into  two  functional  groups. 
First,  the  vortex-lattice  lifting-surface  methodology  is  developed.  An  original  lifting- 
surface  program,  PUF-14,  has  been  written  to  support  the  new  methodology.  This 
first  functional  group  develops  a  method  to  perform  an  analysis  of  the  blade-row 
in  a  specified  inflow.  Areas  of  advancement  include:  wake-adaptive  modeling  in  an 
unsteady  method,  incorporation  of  hub  and  duct  images  in  an  unsteady  method, 
and  Glauret  spaced  panels  which  led  to  a  very  promising  "de-singularized"  Kutta 
condition. 

Second,  the  lifting-surface  method  is  coupled  with  a  RANS  code.  Our  ONR  spon- 
sor has  provided  a  three-dimensional,  steady  RANS  code.  The  RANS  code  is  used 
to  obtain  the  inflow  velocity  field,  which  the  lifting-surface  code  uses  to  calculate  the 
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Figure  1-3:  Typical  maneuvering  transient.   Courtesy  of  SNAME. 

unsteady  forces.  Unsteady  forces  are  generated  due  to  rotating  the  propeller  in  a 
spatially-varying  inflow.  Time- averaged,  but  spatially-varying  body  forces  are  intro- 
duced into  a  three-dimensional  volume  to  represent  propulsor  stages  in  the  RANS 
flow  field.  The  entire  RANS  flow  field  responds  to  the  blade-row  presence.  In  turn, 
the  RANS  flow  field  is  used  again  for  the  lifting-surface  analysis  of  the  blade-row.  By 
alternately  updating  the  lifting-surface  and  the  RANS  solutions,  the  blade-row  forces 
and  RANS  flow  field  converge  to  the  appropriate  solution.  Areas  of  advancement 
include:  determination  of  the  proper  velocities  to  obtain  the  correct  volumetric  effec- 
tive inflow,  determination  of  the  proper  velocities  to  obtain  the  correct  time-average 
forces  and  detailed  study  of  the  relationship  of  time-average  velocities  and  forces  in 
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Figure  1-4:  Typical  steady  inflow  during  maneuvering  transient  is  at  an  angle  (5  with  respect  to  the 
body  centerline. 


Figure  1-5:  Typical  steady  inflow  during  depth  transient  is  at  an  angle  a  with  respect  to  the  body 
centerline. 


both  the  RANS  and  the  potential-flow  domains. 

The  new  treatment  of  unsteady  force  calculations  should  greatly  improve  propulsor 
prediction  capabilities.  In  practice,  past  and  current-day  efforts  treat  unsteady  forces 
with  simplified  assumptions  such  as  cylindrical  propulsor  geometry  and  extremely 
limited  body/propulsor  interaction,  if  any  at  all.  Competing  with  the  new  coupled 
approach  are  fully  three-dimensional,  unsteady  RANS  formulations,  which  are  being 
developed  at  other  institutions.  While  unsteady  RANS  offers  great  potential,  the 
computing  burden  is  enormous  and  not  yet  practical  in  modern  applications.  The  new 
treatment  developed  in  this  thesis  is  believed  to  be  practical  in  both  computational 
load  and  in  representing  the  physical  hydrodynamic  characteristics  of  today's  complex 
propulsors. 
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Chapter  2 

Formulation  of  the  Lifting-Surface 
Problem 

2.1  Fundamental  Assumptions 

The  propeller  is  assumed  to  be  a  set  of  thin  blades  arranged  symmetrically  about 
a  common  axis.  No  restriction  is  made  on  the  blade  shape.  The  propeller  rotates 
with  constant  angular  velocity  in  an  unbounded  fluid.  Hub  and  duct  effects  are 
represented  with  blade  images.  The  onset  flow  is  permitted  to  be  a  function  of  all 
spatial  coordinates. 

The  fluid  is  assumed  to  be  incompressible  and  the  flow  field  irrotational1  except 
on  the  blade  and  in  the  trailing  vortex  wake  sheets.  The  blade  boundary  layer  and 
shed  vortex  wake  thickness  is  assumed  to  be  thin  so  that  the  fluid  rotation  due  to  the 
propeller  is  confined  to  a  thin  layer.  The  fluid  is  treated  as  inviscid  except  for  some 
empirical  corrections  associated  with  the  propeller  blade  drag. 

2.2  Boundary  Value  Problem 

In  the  fluid  domain,  the  velocity,  V,  must  satisfy 

V-V  =  0  (2.1) 


To  be  precise,  the  prescribed  inflow  field  may  be  rotational,  but  vortical  interaction  with  the  potential  flow  field 
induced  by  the  blade  singularities  is  not  treated  explicitly  in  the  blade  solution,  but  rather,  is  accounted  for  by 
interaction  with  a  RANS  code. 
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or.  with  the  assumption  of  the  previous  section, 

V2<?  =  0  (2.2) 

where  <■>  is  the  velocity  potential  defined  by  V  —  V0.  In  a  propeller-fixed  reference 
frame, 

V-n  =  0  (2.3) 

on  the  camber  surface.  At  large  distances  from  the  propeller 

V  ->  U  (2.4) 

where  U  is  the  specified  onset  flow. 

Additionally.  Kelvin's  theorem  for  the  conservation  of  circulation  is  explicitly  in- 
voked on  the  blade  and  in  the  wake  sheet.  The  Kutta  condition  is  implicitly  satisfied 
on  the  trailing  edge  of  each  blade.  Further,  since  a  wake  sheet  is  unable  to  support 
a  pressure  jump  across  it.  the  vorticity  in  the  wake  must  align  with  the  velocities  in 
the  wake  to  remain  force  free. 

2.3     Singularity  Distribution 

The  solution  of  the  boundary  value  problem  is  approached  by  distributing  singularities 
on  the  mean  camber  surface  of  the  propeller  and  on  the  vortex  wake.  Thus,  the  field 
equation,  V24>  =  0,  is  immediately  satisfied. 

A  distribution  of  sources  on  the  camber  surface  is  used  to  generated  the  jump  in 
normal  velocity.  The  source  strengths  are  obtained  by  stripwise  application  of  thin 
wing  theory.  The  jump  in  tangential  velocity  across  the  camber  surface  is  provided 
by  a  distribution  of  vorticity  on  the  camber  surface  and  in  the  wake  sheet.  The 
vorticity  strengths  are  obtained  by  imposition  of  the  boundary  conditions  stated 
above.  This  leads  to  a  singular  integral  equation  over  the  blade  and  wake  surface.  By 
the  numerical  scheme  described  later,  the  boundary  problem  is  reduced  into  a  system 
of  linear  algebraic  equations. 
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Chapter  3 

Mathematical  Modeling  of  the 
Propeller  System 

3.1      Blade  Geometry 

Propeller  blades  are  traditionally  defined  by  tabular  data  consisting  of  radial  distribu- 
tions of  geometrical  quantities  such  as  pitch,  rake  and  skew,  and  by  chordwise  distri- 
butions of  camber  and  thickness  defined  along  circular  cylinders.  However,  problems 
arise  due  to  interpolation  and  smoothing  during  each  step  in  the  design  and  man- 
ufacturing process.  This  problem  is  intensified  as  the  design  process  becomes  more 
complex,  where  geometry  must  be  passed  among  a  large  number  of  hydrodynamic 
and  structural  analysis  codes. 

To  uniquely  define  the  blade  geometry,  B-spline  surfaces  are  used.  A  B-spline 
surface  can  be  interrogated  at  an  arbitrarily  fine  mesh  of  points  to  define  the  blade. 
Therefore,  designers  and  manufacturers  each  can  use  the  same  B-spline  to  define  the 
geometry  to  the  required  accuracy  of  their  respective  calculations  [35]. 

The  propeller  blade  surface  is  described  by  means  of  a  control  polygon  net  of  B- 
spline  vertices.  The  number  of  control  polygon  vertices  required  to  define  the  camber 
surface  of  a  propeller  blade  can  be  quite  small,  with  a  7  x  7  grid  being  satisfactory 
in  many  applications.  The  B-spline  surface  is  easy  to  manipulate  during  the  design 
process,  such  as  when  designing  with  PBD-14  [39],  because  each  B-spline  vertex  has 
a  local  effect  as  shown  in  figure  3-1. 
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Surface  Grid 


Displaced  Vertex 


Control  Net  and  Surface  Grid 


Figure  3-1:  The  effect  on  blade  shape  of  moving  one  B-spline  control  vertex. 

Pre-  and  post-processors  enable  the  use  of  the  designer's  description  of  the  blade 
surface,  i.e  pitch,  rake,  skew.  The  designer's  descriptions  are  necessary  when  commu- 
nicating with  older  codes/tools.  The  tools  used  for  these  conversions  are  described  in 
the  Propeller  Blade  Design  (PBD-14)  Manual  [39]. 

The  blade  coordinate  system  used  in  this  work  closely  follows  that  described  in 
[29]  and  [14].  A  Cartesian  coordinate  system  is  fixed  to  the  propeller  with  the  x- 
axis  coaxial  with  the  propeller  hub  and  pointing  in  the  streamwise  direction.  The 
y-axis  may  be  oriented  at  any  convenient  angle  to  the  propeller.  The  corresponding 
cylindrical  coordinates  are  defined  with  r2  =  y2  +  z2  and  the  blade  rotation  angle,  6, 
is  measured  from  the  y-axis  in  a  right-handed  sense.  Figure  3-2  illustrates  the  blade 
geometry  and  the  coordinate  system. 

3.2     Discretization  of  Blade  Singularity  Distribution 

The  method  of  singularity  distribution  is  one  of  the  most  powerful  techniques  for 
the  solution  of  the  fluid  flow  problem.  The  boundary  value  problem  formulated  in 
chapter  2  can  be  divided  into  two  distinct  groups:  the  source  singularities  and  the 
vortex  singularities. 
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Downstream 


(Rotates 
with  Blade) 


Figure  3-2:  Pictorial  of  the  blade  coordinate  system.  Note  that  a  right-handed  propeller  rotates  in 
a  direction  of  negative  angle. 


The  continuous  distribution  of  sources  and  vortices  is  replaced  by  an  array  of 
M  x  A  concentrated  straight-line  elements  of  constant  strength.  The  end  points 
of  the  elements  are  located  on  the  camber  surface.  The  exact  arrangement  of  the 
points  has  been  the  subject  of  much  experimentation.  The  chosen  spacing  on  the 
camber  surface  follows  the  method  elected  by  Greeley  and  Kerwin  [14].  Figure  3-3 
illustrates  the  straight-line  representation  of  the  continuous  singularity  distributions 
for  a  M  x  Ar  =  11  x  5  blade  grid. 

3.2.1      Source  Distribution 

The  source  singularities  are  distributed  to  represent  the  jump  in  normal  velocity  across 
the  camber  surface.  At  the  outset,  we  assume  the  source  strength  distribution  is 
independent  of  time,  and  that  its  spatial  distribution  may  be  derived  from  a  stripwise 
application  of  thin-wing  theory  at  each  radius.  These  assumptions  are  consistent  with 
Kerwin  and  Lee  [29]. 
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Figure  3-3:  This  figure  shows  the  blade  grid  and  wake  grid  points  connected  by  a  mesh.  The 
blade  control  points  are  also  shown.  The  boundary  value  problem  solves  for  the  blade  circulation 
distribution  which  simultaneously  nulls  the  normal  velocity  at  every  control  point  on  every  blade 
surface. 


Strictly  speaking,  the  source  strengths  are  dependent  on  the  local  inflow  and  will 
therefore  vary  in  time.  However,  by  using  the  mean  inflow  to  set  the  source  strengths. 
the  variations  with  time  from  the  mean  value  will  be  small.  Additionally,  the  blade 
thickness  contribution  to  the  boundary  value  problem  is  secondary  considering  both 
the  mean  and  fluctuating  blade  loading.  In  recognition  of  the  secondary  influence 
of  blade  thickness,  it  is  a  reasonable  assumption  that  the  source  strengths  are  time 
invariant.  Thus,  the  source  strength  and  their  contribution  to  the  boundary  value 
problem  are  solved  only  once. 

3.2.2     Vortex  Distribution 

The  jump  in  tangential  velocity  across  the  camber  surface  is  provided  by  a  distribution 
of  vorticity  on  the  camber  surface  and  in  the  wake  sheet.  As  shown  in  figure  3-3,  the 
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vortex  strength  is  a  vector  lying  on  the  surface  and  may  be  resolved  into  components 
along  two  arbitrarily  assigned  directions  on  the  surface.  The  vortex  distribution  on 
the  blade  is  chosen  to  be  resolved  into  "spanwise"  and  "chordwise"  components  while 
the  corresponding  components  in  the  wake  are  termed  '"shed"  and  "trailing"  vorticity. 
In  general,  the  shed  vorticity  is  a  result  of  temporal  variations  of  the  spanwise  vorticity 
strength  on  the  blade. 

An  implicit  Kutta  condition  is  established  by  arranging  the  blade  grid  and  the 
control  points  such  that  a  control  point  lies  on  the  blade  trailing  edge.  Thus,  when 
the  condition  of  V  -ft  =  0  is  satisfied  exactly  the  blade  trailing  edge,  the  Kutta 
condition  is  simultaneously  satisfied.  The  implicit  Kutta  condition  has  been  shown 
to  work  in  Greeley  and  Kerwin  [14]  and  Breslin.  et  al.  [4]. 

3.3      Geometry  of  Transition  and  Ultimate  Wakes 

The  geometry  of  the  trailing  vortex  wake  has  an  important  influence  on  the  accuracy 
of  the  calculation  of  induced  velocities  on  the  blade.  The  wake  alignment  procedure 
employed  by  Greeley  and  Kerwin  [14]  convects  the  trailing  vortices  along  axisymmet- 
ric  surfaces  constructed  from  a  user-supplied  tip  contraction  angle  and  ultimate  wake 
radius.  The  present  scheme  follows  that  used  in  PBD-14  and  allows  the  wake  to  con- 
form to  the  actual  inflow  velocity  field.  Thus,  the  wake  conforms  to  the  body  shape 
and  interior  duct  surface,  if  present,  and  follows  the  circumferential-mean  stream- 
tubes. 

The  wake  is  divided  into  two  regions: 

1.  Transition  Wake.  The  transition  wake  region,  wherein  all  contraction  and 
roll  up  of  the  trailing  vortex  sheet  is  assumed  to  occur,  is  built  upon  the  actual 
inflow  velocity  field.  The  transition  wake  is  an  extension  of  the  blade  vortex- 
lattice  grid.  The  axial  extent  of  the  transition  wake  is  specified  by  the  user. 
The  shed  vorticity  in  the  transition  wake  is  assumed  to  decay  in  strength  as  it 
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is  washed  downstream. 

2.  Ultimate  Wake.  The  ultimate  wake  comprises  the  region  from  the  end  of  the 
transition  wake  to  a  point  infinitely  far  downstream.  The  trailing  vortex  lines 
from  each  blade  are  merged  into  an  infinite-bladed  helical  vortex  whose  radii 
match  that  of  the  downstream  end  of  the  trailing  vortex  lines.  The  total  cir- 
culation of  the  helical  vortex  is  determined  by  the  average  of  the  shed  vorticity 
at  each  individual  trailing  vortex.  Efficient  closed-form  expressions  for  the  ve- 
locity induced  by  infinite-bladed  helical  vortices  were  developed  by  Hough  and 
Orel  way  [15].  Leibman  [32]  demonstrated  that  the  infinite-bladed  approximation 
introduced  negligible  errors  provided  that  the  transition  wake  length  was  at  least 
one  propeller  radius. 

The  transition  wake  grid  is  implicitly  specified  by  the  grid  density  on  the  blade 
and  by  the  time  stepping  variables.  The  trailing  vortices  in  the  transition  wake  are 
extensions  of  the  chordwise  vortices  on  the  blade.  The  discretized  representation  of 
the  vortex  sheet  then  consists  of  M  concentrated  trailing  vortex  lines  whose  trailing 
edge  coordinates  match  the  corresponding  values  of  the  chordwise  vortices  on  the 
blade. 

The  transition  wake  region  contains  a  set  of  Nj  x  M  shed  vortex  segments,  where 
Nj  is  the  number  of  shed  vortex  segments  in  the  transition  wake.  The  segments 
occupy  the  region  from  the  blade  trailing  edge  to  the  start  of  the  ultimate  wake 
The  total  inflow  which  convects  the  transition  wake  elements  can  vary  with  angular 
position  as  seen  in  figure  3-4.  However,  as  a  simplifying  assumption,  the  wake  is  built 
upon  the  circumferential-mean  inflow.  The  wake  geometry  has  no  dependence  on  the 
blade  angular  position;  thus,  the  wake  geometry  remains  unchanged  relative  to  the 
blade.  Thereby,  a  single  wake  geometry  can  be  used  to  represent  all  transition  wakes 
at  all  angular  positions.  While  this  is  artificial,  it  is  deemed  necessary  considering 
the  scope  of  the  project  and  is  expected  to  add  a  negligible  error  in  the  calculations. 
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Figure  3-4:  This  figure  shows  a  center-body  representation  and  the  blade  grid  overlayed  with  t he 
input  axial  velocity. 


Tracking  of  individual  wakes  trajectories  could  be  added  in  a  future  revision. 

The  angular  increment,  59,  and  the  time  increment,  5t,  are  related  through  the 
equation 

50  =  u)6t  (3.1) 

where  lj  is  the  rotation  rate  of  the  propeller. 

Given  the  non-dimensional  assumptions  of  Vs  =  1  and  Rbiadetip  —  15  the  convection 
parameter  which  governs  the  time  interval  between  each  time  step  is 

2JS 


St  = 


Ne 


(3.2) 


where  Js  is  the  advance  coefficient  and  Ng  is  the  number  of  time  steps  per  revolution. 
Thus,  the  convection  of  the  shed  vorticity  is  related  to  all  components  of  the  inflow 
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velocity  and  to  the  time  increment  implicitly  specified  by  the  input  parameters.  This 
improves  upon  the  linear  convection  used  in  Kerwin  and  Lee  [29]  and  Keenan  [24]. 
[23].  Additionally,  by  design,  there  exists  an  explicit  correspondence  of  the  shed 
vortex  geometry  to  the  time  step  interval,  the  nth  shed  vortex  at  time  /  will  have  the 
same  strength  as  the  (n  —  l)th  vortex  at  time  (t  —  1). 

While  the  shed  vorticity  is  convecting  downstream,  it  is  hypothesized  that  the 
shed  vorticity  becomes  highly  disorganized  within  a  short  distance  behind  the  trailing 
edge  [29].  In  order  to  simulate  this  dissipation  of  shed  vorticity,  a  decay  factor  7(P) 
is  applied  to  the  shed  vorticity  strength  given  by 

7(P)  =  2P3-3P2  +  1  0<P<1  (3.3) 

where  P  is  the  fractional  distance  along  the  vortex  sheet  from  the  blade  trailing  edge 
to  the  ultimate  wake  [23].  In  this  way,  the  shed  vortex  strength  is  made  to  approach 
the  average  value  for  that  particular  M  set  of  vortices.  Conveniently,  the  ultimate 
wake  approximation  uses  this  average  value  to  estimate  the  strength  of  the  ultimate 
wake.  Thereby,  a  smooth  transition  of  vortex  strength  is  maintained  from  the  blade 
through  the  transition  wake  and  into  the  ultimate  wake. 

3.4     Representation  of  Propeller  and  Wake  Vorticity 

The  propeller  is  assumed  to  have  several  blades  which  are  identical  in  shape  and 
evenly  spaced  around  a  common  axis.  Some  earlier  analysis  techniques  solved  the 
boundary  value  problem  by  solving  the  problem  on  only  one  blade.  1  However,  for 
this  solution  method,  the  boundary  value  problem  includes  all  blades.  Thus,  every 
control  point  on  every  blade  is  simultaneously  solved. 

The  effects  of  the  vortex  elements  are  taken  into  account  by  conceptualizing  the 
vorticity  as  closed  loops  which  automatically  satisfy  Kelvin's  theorem.    Figure  3-5 

1  Specifically,  PBD-14  is  an  steady,  axisymmetric  solver  which  solves  one  blade  BVP.  The  influence  of  all  blades  are 
correctly  accounted  for  through  the  influence  function.  Conversely,  PUF-2  used  a  coarse  spacing  on  all  blades  except 
the  key  blade  and  iteratively  approached  a  solution  at  each  time  step. 
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shows  the  conceptual  idealization.  The  influence  of  a  blade  vortex  element  on  the 
velocity  at  a  control  point  is  captured  with  a  Loop  Influence  Function.  The  induced 
velocities  due  to  the  LIF  are  unknown  for  the  current  time  step.  The  LIF  is  placed 
on  the  left-hand-side  of  equation  3.4. 

The  influence  of  a  wake  vortex  element  is  captured  using  a  Wake  Influence  Func- 
tion, WIF .  The  wake  influence  depends  on  the  time  history  of  blade  loading,  thus  the 
induced  velocities  on  a  control  point  are  known  quantities.  Consequently,  the  wake 
vorticity-induced  velocities  are  placed  on  the  right  hand  side  of  equation  3.4. 

The  LIFs  and  WIFs  are  developed  as  the  velocity  induced  on  a  control  point  due 
to  a  unit  strength  vortex  loop.  Lee  [31]  and  Keenan  [24]  show  the  construction  of  the 
simultaneous  equations  for  closed  loops.  The  present  work  follows  Lee  and  Keenan 
in  the  derivation  of  loops. 

Unlike  Lee  and  Keenan.  the  present  work  solves  the  boundary  value  problem  on 
all  blades  simultaneously.  The  formulation  of  the  simultaneous  equations,  while  con- 
ceptually similar,  yields  many  more  equations.  An  additional  difference  is  that  the 
present  formulation  uses  cosine  spacing  across  the  chord  of  the  blade  with  an  implicit 
Kutta  condition  at  the  trailing  edge. 

As  expected,  the  LIF  matrix  exhibits  the  structure  for  the  propeller  system.  Equa- 
tions 3.4  and  3.5  show  the  structure  that  exists  in  the  simultaneous  equations.  These 
equations  show  the  influence  functions  after  they  have  been  converted  to  represent 
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Figure  3-5:  The  boundary  value  problem  is  solved  by  conceptualizing  the  vorticity  as  closed  loops 
on  the  blade  and  in  the  wake.  The  ultimate  wake  is  a  "closed"  loop  extending  to  infinity. 


the  normal  component  of  velocity  due  to  a  unit-strength  vortex  loop. 
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(3.4) 


Accounting  for  symmetry  such  as  Blade-on-Blade  which  is  the  same  for  all  blades 
and  Blade-on-Prior-Blade  which  is  the  same  for  all  Blade-on-Prior-Blade  relations. 
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then  the  structure  is  more  evident.  Equation  3.5  assumes  a  five-bladed  propeller  and 
replaces  LI F"  with  .4,  LIFlf  with  B,  etc..  to  use  the  Blade-on--  •  •  relations  which 
makes  the  structure  much  more  clear. 
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3.5     Stepwise  Solution  in  the  Time  Domain 


[r,]  =  -iv -h] 


[3.5) 


With  the  knowledge  of  how  the  propeller  system  is  to  be  modeled,  the  sequence  of 
program  execution  must  be  considered.  Figure  3-6  shows  a  functional  block  diagram 
of  how  the  program  execution  should  proceed.  Note  that  the  simultaneous  equations 
must  be  solved  at  each  time  step.  All  information  that  is  not  time  dependent  is 
calculated  prior  to  entering  into  the  time  domain  calculations. 
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Figure  3-6:  Functional  block  diagram  of  the  PUF-14  methodology 
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Chapter  4 

Determination  of  Blade  Forces 


The  force  and  moment  acting  on  the  propeller  blade  can  be  obtained  by  integrating 
the  pressure  jump  over  the  blade  camber  surface. 

F  =  f  ApndA  (4.1) 

M  =   f  Ap(r  x  n)dA  (4.2) 

where  A/>  is  the  pressure  jump  across  the  camber  surface,  n  is  the  normal  vector 
on  the  camber  surface,  and  r  is  the  position  vector  from  the  origin  to  the  point  of 
integration. 

The  forces  acting  on  a  propeller  blade  can  be  divided  into  four  portions:  the 
pressure  forces  acting  normal  to  the  blade  surface,  the  viscous  forces  acting  tangential 
to  the  blade  surface,  the  leading  edge  suction  force  and  the  force  proportional  to  the 
time  rate  of  change  of  potential,  which  follows  from  the  unsteady  term  in  Bernoulli's 
equation.  The  four  components  of  blade  force  are  discussed  below. 

4.1      Pressure  Forces 

The  force  acting  normal  to  the  blade  surface,  commonly  referred  to  as  the  pressure 
force,  is  obtained  from  JoukowskTs  law.  At  the  control  point  in  a  given  panel  there 
is  a  vortex  density  7  corresponding  to  the  local  pressure  jump  across  the  blade.  The 
Joukowski  force  on  the  panel  is 

dFj  =pdA(V  x7),  (4.3) 
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where  V  is  the  mean  velocity  acting  at  the  control  point  of  the  panel.  d.\  is  the  area 
of  the  panel  and  p  is  the  fluid  density.  The  vortex  density  7  is  related  to  the  jump  in 
tangential  velocity,  2  dv,  across  the  blade  by 

7  =  n  x  2  dv,  (4.4) 

n  being  the  unit  normal  to  the  blade  surface.  This  can  be  used  to  rewrite  the 
Joukowski  force  as 

dFj  =  pd.\  ((V-2f/v)n-  (V-n)2<fv).  (4.5) 

The  jump  in  tangential  velocity,  2  dv,  is  obtained  by  taking  the  gradient  of  the 
potential  jump  across  the  vortex  sheet 

2dv  =  Vsn.  (4.6) 

Vs  is  the  gradient  operator  on  the  blade  surface  and  //  is  the  potential  jump  across 
the  vortex  sheet.  In  a  vortex-lattice  representation,  the  potential  jump  at  the  control 
point  is  equal  to  the  vortex  loop  strength  around  the  panel,  so  the  jump  velocity  is 
obtained  directly  by  differentiating  the  vortex  loop  strengths  representing  the  loaded 
blade  surface. 

4.2  Viscous  Forces 

Viscous  forces  are  computed  using  an  empirical  sectional  drag  coefficient.  This  coef- 
ficient, Cdv,  is  supplied  by  the  user.  An  element  of  viscous  force  acting  on  one  panel 
of  the  blade  surface  is  then  computed  as 

dFv=l-PV\V\CdvdA.  (4.7) 

The  total  drag  force  is  then  obtained  by  summing  the  elemental  forces  over  all  panels. 

4.3  Leading  Edge  Suction  Forces 

A  blade  which  is  not  acting  at  ideal  angle  of  attack  (shock-free  entry)  has  an  additional 
force  acting  on  it — the  leading  edge  suction — which  is  a  consequence  of  representing 
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the  blade  as  a  thin  wing.  The  suction  force  on  a  segment  of  the  leading  edge  is 
computed  as 

dFL  =  -l-7rpC2sb,  (4.8) 

where  b  is  a  unit  vector,  orthogonal  to  the  leading  edge,  lying  on  the  camber  surface 
of  the  blade  and  Cs  is  the  suction  coefficient.  This  coefficient  is  obtained  from  the 
limit 

C.  =  limv^(7M-t),  (4.9) 

s->0 

where  t  is  a  unit  vector  tangent  to  the  leading  edge,  and  s  is  the  curvilinear  distance 
along  the  camber  surface  from  the  leading  edge.  The  total  suction  force  is  obtained 
by  summing  the  elemental  dF^s  over  the  length  of  the  leading  edge  [25]. 

4.4     Unsteady  Velocity  Potential 

As  the  blade  passes  through  a  region  where  the  inflow  changes,  the  bladeexperiences 
an  added-mass  force.  This  force  is  proportional  to  the  time  derivative  of  the  veloc- 
ity potential  and  follows  from  the  unsteady  term,  equation  4.10,  from  Bernoulli's 
equation. 

d¥A  =  pnj{6+  -  <jT)dA  (4.10) 

First,  recall  that  the  source  strengths  are  assumed  to  the  independent  of  time,  so  that 
only  the  potential  due  to  the  vortices  is  required.  Since  the  jump  in  potential  across 
the  camber  surface  in  the  continuous  case  is 

0+00  -4T(a)=  Is    7(CMC  (4-11) 

Jl.e. 

we  can  replace  4.10  by  the  final  form 

dFA=pn±  f    fiOdCdA  (4.12) 

at  jl.e. 

The  time  derivative  in  4.12  may  be  obtained  by  numerical  differentiation  of  the  dis- 
crete vortex  strengths  obtained  at  several  discrete  time  steps.  Since  the  calculations 
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of  the  forces  can  be  made  after  the  results  of  the  propeller  have  been  time-domain 
solution,  the  future  is  already  known:  thus,  permitting  the  use  of  a  very  accurate 
central  five-point  differentiation  formula. 
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Chapter  5 

Validation  using  a  Specified 

Effective  Inflow  H 

5.1      Convergence 

There  are  many  parameters  that  must  converge  to  yield  the  final  converged  answer. 
The  first  parameter  studied  is  the  convergence  of  the  solution  for  various  blade  grid 
resolutions.  Propeller  4577  grid  resolution  is  varied  from  a  coarse  resolution  of  S  x  S 
lattice  to  a  much  finer  lattice  of  36  x  36.  The  chosen  inflow  is  described  in  section  5.3. 
The  results  are  plotted  in  figure  5-1.  The  figure  shows  good  convergence  of  Ay  and 
I\q  for  a  relatively  coarse  grid.  In  practice,  the  desired  grid  is  one  which  gives  the 
desired  accuracy  with  the  least  computational  effort.  A  grid  of  around  12  x  12  provides 
a  reasonable  compromise. 

The  next  convergence  study  centers  on  unsteady  analysis.  Propeller  41  IS  is  used  to 
test  the  results  of  PUF-14  for  convergence  in  an  unsteady  analysis.  The  chosen  wake  is 
the  4118  wake,  published  in  [29],  with  slight  modifications  at  the  hub  and  tip  to  ensure 
the  wake  extends  to  the  blade  extremes.  The  wake  is  essentially  described  as  one  plane 
of  flow  data  which  varies  circumferentially  and  with  radius.  This  particular  4118  wake 
was  experimentally  made  by  superimposing  wakes  that  have  strong  harmonic  content 
in  the  third  and  fourth  harmonic. 

The  solved  circulation  for  the  entire  time  domain  solution  is  shown  in  figure  5-2. 
When  the  ''starting  vortex"  is  first  generated,  there  is  a  strong  effect  on  the  solved 
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Figure  5-1:  Convergence  of  the  forces  with  increasing  blade  grid  resolution. 

circulation.  As  the  "starting  vortex"  washes  downstream,  its  strong  local  effects  are 
mitigated  which  allows  the  solved  circulation  to  approach  its  converged  value  for  each 
particular  angular  position. 

The  next  significant  parameter  to  test  for  convergence  is  the  size  of  the  time  step. 
This  parameter,  called  the  number  of  time  steps  per  revolution  or  Ne,  controls  the 
size  of  the  angular  rotation  per  time  step  which,  in  turn,  controls  the  downstream 
convection  distances  of  each  vortex  loop.  Figure  5-3  shows  the  error  associated  with 
various  time  step  sizes.  The  plotted  parameter  is  the  circulation  at  the  0.7  r/R  of  the 
blade.  The  circulation  is  strongly  influenced  by  the  harmonics  in  the  inflow  velocity 
field.  The  converged  solution  for  circulation  using  210  time  steps  per  revolution  is 
taken  as  the  reference  circulation.  Thus,  the  plotted  values  represent  the  percent 
error  associated  with  a  given  Ng. 

Note  that  in  figure  5-3  the  average  magnitude  of  error  does  not  improve  signifi- 
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Figure  5-2:  The  solved  circulation  is  shown  for  the  entire  time  domain  solution  calculated  by  PUF-14. 
Note  that  only  the  converged  solution  correctly  represents  the  boundary  value  problem. 


cantly  above  60  N$.  Thus,  a  value  of  60  N$  is  a  good  compromise  between  speed 
and  accuracy.  Figure  5-4  shows  a  direct  comparison  of  two  values  of  the  Nq.  60  and 
210.  These  two  cases  show  that  the  solved  blade  circulation  tracks  closely  for  the  two 
values  of  Ne- 

In  viewing  figure  5-3,  there  is  obviously  a  periodicity  to  the  error  associated  with 
the  convergence  of  the  time  step  size.  The  maximum  error  tends  to  occur  where  the 
solved  circulation  has  the  greatest  slope.  Therefore,  a  distinction  in  the  error  associ- 
ated with  amplitude  and  phase  may  be  seen  through  a  comparison  of  the  harmonic 
analysis. 

Figures  5-5  and  5-6  show  the  decomposition  of  the  convergence  error  associated 
with  the  force  calculations  into  amplitude  and  phase.  The  amplitude  is  fairly  well 
converged.    However,  the  phase  is  not  so  settled  and,  in  fact,  becomes  quite  poor 
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Figure  5-3:  Convergence  error  associated  with  the  angular  size  of  the  time  step.  Arg 
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Figure  5-4:  Comparison  of  two  angular  sizes  of  the  time  step,  Ng 

in  the  higher  harmonics.  The  errors  in  phase  are  somewhat  expected  due  to  the 
arbitrary  choices  in  converting  the  continuous  boundary  value  problem  into  a  discrete 
problem  of  simultaneous  equations.  The  arbitrary  choices  become  less  important  as 
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the  discrete  time  step  becomes  smaller  to  approach  the  continuous  case.  Therefore, 
when  accurate  phase  information  is  need  at  the  higher  harmonics,  smaller  time  steps 
are  needed. 
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Figure  5-5:  Comparison  of  harmonic  phases  of  various  time  step  sizes.  Xg 
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Figure  5-6:  Comparison  of  harmonic  amplitudes  of  various  time  step  sizes,  Ng 
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5.2  Validation  Against  Other  Numerical  Methods 

The  stand-alone  mode  of  PUF-14  was  validated  in  an  effective  wake  by  comparing  to 
the  results  of  PUF-2.1  and  PBD-14.2.8  in  the  identical  wakes,  respectively.  Since  these 
two  codes,  PUF-2.1  and  PBD-14.2.8,  have  been  in  existence  for  many  years,  they  have 
been  validated  both  at  the  initial  code  validation  phase  and,  later,  validated  against 
real-life  designs.  Thus,  the  PUF-14  computer  program  can  be  validated  against  these 
two  pre-existing  programs  for  the  specialized  cases  discussed  below. 

5.3  Steady  Analysis 

PBD-14.2.8  has  been  validated  in  such  references  as  [11],  [2],  [1].  Comparison  against 
PBD-14  will  be  used  to  validate  PUF-14  in  the  steady  mode  of  analysis.  The  chosen 
test  case  uses  the  propeller  4577  from  Huang  [IS].  The  effective  inflow  is  approximated 
by  using  the  total  inflow  from  a  converged  axisymmetric  RANS  flow  field  for  the 
Huang  body  one.  While  the  effective  inflow  may  not  be  correct  for  this  propeller,  it 
provides  a  rough  estimate  of  an  effective  flow  field  and  allows  a  direct  comparison  of 
the  results  of  the  two  codes. 

5.3.1  Solved  Blade  Circulation 

The  two  cases  shown  in  figure  5-7  illustrate  the  agreement  between  PUF-14  and  PBD- 
14.2.8.  Additionally,  cases  run  using  various  combinations  of  parameters,  including 
stator  blade-rows,  all  show  nearly  exact  agreement  between  PBD-14. 2. S  and  PUF-14 
for  steady  axisymmetric  test  cases. 

5.3.2  Forces 

The  two  cases  shown  in  figure  5-7  have  a  steady  value  of  thrust  and  torque.  These 
values  are  shown  in  tables  5.1  and  5.2.  The  results  are  exact,  indicating  that  PUF-14 
is  correctly  solving  the  special  case  of  steady,  axisymmetric  propeller  analysis.  As 
stated  above,  additional  cases  all  show  similar  agreement. 
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Figure  5-7:  Comparison  of  the  Non-Dimensional  Circulation  versus  Radius 


A'j 

Kq 

PBD-14.2.8 
PUF-14 

0.0838 
0.0838 

0.1920 
0.1920 

Table  5.1:  Steady  axisymmetric  force  comparison  for  the  case  with  thickness,  hub,  duct.  6  blades 
with  a  10x10  grid 


Kt 

i<q 

PBD-14.2.8 
PUF-14 

0.0567 
0.0567 

0.1219 
0.1219 

Table  5.2:  Steady  axisymmetric  force  comparison  for  the  case  with  no  hub,  duct  nor  thickness,  3 
blades  with  a  12x12  grid 


5.4     Unsteady  Analysis 

PUF-2.1  has  been  validated  by  many  authors  [31],  [29],  [26].  Thus,  comparisons 
with  PUF-14  will  suffices  for  validation  of  correctly  solving  the  specialized  problem 
of  circumferential-  and  radial-varying  inflow  but  no  variation  along  the  propeller  axis 
and  no  hub  nor  duct.  The  chosen  test  case  is  propeller  4118  as  discussed  earlier  in 
this  chapter. 
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The  most  significant  difference  between  the  codes  leading  to  the  different  solved 
circulation  is  the  lattice  arrangement  on  the  blades  and  in  the  wake.  Figure  5-S  shows 
the  key  blade  and  the  transition  wake  for  the  two  codes.  Recall  that  PUF-14  uses 
identical  blade  lattices  and  identical  transition  wakes  for  all  blades.  PUF-2.1  uses  a 
much  coarser  lattice  on  blades  other  than  the  key  blade. 


Figure  5-8:  The  leftmost  figure  shows  PUF-2.1  key  blade  and  transition  wake  lattice  arrangement. 
The  rightmost  figure  shows  PUF-14  key  blade  and  transition  wake  lattice  arrangement.  Differences 
in  lattices  lead  to  differences  in  solved  blade  circulation  shown  later. 


More  specifically,  PUF-2.1  uses  constant  spacing  in  the  chordwise  direction  while 
PUF-14  uses  cosine  spacing.  Therefore,  PUF-14  resolves  more  of  the  leading  and 
trailing  edge  gradients.  Another  lattice  difference  is  in  the  wake  models:  PUF-2.1 
only  propagates  about  a  quarter  of  a  propeller  diameter  aft  while  PUF-14  propagates  a 
user-specified  distance  aft.  Another  difference  exists  in  the  specification  of  the  general 
shape  for  the  transition  wake  and  in  the  ultimate  wakes.  Still  another  difference 
exists  in  the  specification  of  the  blade  geometry  and  lattice.  That  is,  PUF-2.1  uses 
cylindrical  geometry  parameters  to  built  its  lattice  while  PUF-14  uses  a  B-spline 
geometry  description  overlayed  with  a  vortex  lattice  built  on  the  circumferential- 
mean  streamtubes. 
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5.4.1      Solved  Blade  Circulation 

Figure  5-9  shows  the  agreement  of  the  blade  circulation  at  the  0.7  r/R  position.  Some 
differences  of  circulation  are  expected  due  to  the  inherent  differences  in  the  solution 
methodology  between  the  two  codes. 


-0 PUF14TIMESOLN 

— PUF2TIMESOLN 


100  200  300 

Key  Blade  Position  (degrees) 


Figure  5-9:  Comparison  of  the  Non-Dimensional  Circulation  at  0.7  r/R 

5.4.2      Forces 

The  cases  shown  in  figures  5-10  and  5-11  compare  the  rr,  y,  and  z  forces  and  mo- 
ments, respectively,  on  a  blade  as  it  rotates.  Note  that  the  xyz  coordinate  system 
rotates  with  the  blade.  As  discussed  above,  exact  agreement  of  the  solved  circulation 
and  forces  is  not  expected.  In  the  figures,  there  is  evidence  of  a  small  error  in  the 
magnitude  of  forces/moments  and  a  small  error  in  the  phase  of  the  forces/moments. 
However,  there  is  significant  agreement  between  the  two  solution  methods  which  leads 
to  the  conclusion  that  PUF-14  is  correctly  solving  the  specialized,  unsteady  propeller 


47 


boundary  value  problem  of  cylindrical  flow. 

5.5     Validation  Against  Experiment 

5.5.1  Experimental  Setup 

Experiments  are  described  in  reference  [21]  in  which  the  propeller  4679  is  used  to 
obtain  unsteady  pressure  distributions  in  inclined  flow.  The  propeller  was  designed 
to  model  a  controllable-pitch  (CP)  propeller  used  for  high-powered  ships.  A  picture 
of  propeller  4679  is  shown  in  figure  5-12. 

The  experiment  was  conducted  on  the  DTIVSRDC  Carriage  V.  The  propeller  was 
driven  from  downstream  while  the  carriage  supporting  the  propeller  advanced  into 
the  flow.  The  comparisons  chosen  for  validation  of  PUF-14  use  the  test  in  which  the 
propeller  shaft  was  inclined  at  7.5  degrees  downward  from  the  direction  of  the  carriage 
advance  so  that  the  propeller  operated  in  inclined  flow.  Two  advance  coefficients  have 
been  used  for  comparison  of  PUF-14  and  this  experiment. 

No  flow  field  measurements  were  given  so  the  geometric  flow  field  is  used  to  es- 
timate the  propeller  inflow.  The  inclined  flow  is  represented  by  once-per- revolution 
variation  in  the  radial  and  tangential  flow  components.  As  a  further  refinement  of  the 
flow  field,  the  prescribed  geometrical  inflow  was  modified  with  a  crude  adjustment 
for  the  presence  of  the  propeller  shaft.  An  axial  potential  flow  approximation  was 
included  by  modeling  the  propeller  shaft  as  a  Rankine  ovoid  centered  on  the  propeller 
plane.  A  similar  correction  was  applied  for  the  cross-flow  directions  using  the  poten- 
tial flow  around  a  two-dimensional  cylinder.  Convection  velocities  were  used  for  wake 
alignment  as  computed  from  PSF-2  [14]. 

5.5.2  Comparisons  with  Experiment 

Figures  5-13  and  5-14  show  the  mean  pressure  distribution  on  one  blade  of  DTMB 
4679  at  r/R  =  0.7  for  Js  of  1.078  (design  Js)  and  0.719,  respectively.  The  figures 
compare  PUF-14,  PUF-10.3  and  the  experimental  results.    PUF-10  is  an  unsteady, 
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low-order  perturbation-potential-based  panel  method  incorporating  an  unsteady  it- 
erative pressure  Kutta  condition  at  the  blade  trailing  edge.  The  trailing  vortex  wake 
has  a  frozen  geometry  based  upon  user  supplied  inflow  inputs  and  convection  veloc- 
ities. PUF-10  was  given  the  same  specified  inputs  as  PUF-14.  PUF-10  provides  a 
check  of  the  trends  expected  for  the  pressure  distribution.  The  experimental  results 
provide  point  values  expected  from  the  calculations.  PUF-10. 3  and  PUF-14  follow 
the  same  basis  shape  while  nearly  matching  the  experimental  data  points. 

Figures  5-15  through  5-1S  show  a  comparison  of  the  first  harmonic  amplitudes 
and  phases.  The  figures  compare  the  experimental  results  with  the  calculational 
results  of  PUF-14.  The  PUF-10  results  are  not  in  agreement  with  either  PUF-14 
or  the  experiment  so  those  results  are  not  shown  here.1  A  definite  agreement  is 
established  between  PUF-14  calculations  and  the  experiment.  The  two  agree  for  on- 
and  off-design  advance  coefficients  in  both  the  first  harmonic  amplitude  and  the  first 
harmonic  phase. 

The  reasons  for  the  differences  that  do  exist  are  unclear.  One  reason  is  likely  due 
to  the  source  strength  representation  of  blade  thickness.  In  PUF-14.  the  line  sources 
used  to  represent  the  blade  thickness  are  solved  at  the  onset  of  the  solution  using  the 
circumferential-mean  inflow.  Therefore,  the  source  strengths  are  invariant  in  time  and 
space.  This  assumption  means  that  the  once-per-revolution  variation  in  tangential 
velocity  on  the  blade,  which  linearly  determines  the  source  strengths,  is  not  captured. 
Thus,  the  net  impact  is  that  velocities  variations  at  higher  harmonics  are  not  totally 
correct.  That  is,  the  velocity  does  not  include  the  perturbation  velocities  which  would 
be  due  to  the  source  strength  variation  as  the  blade  rotates.  This  simplification 
does  not  significantly  affect  the  force  calculations  since  forces  are  derived  from  a 
pressure  difference.  However,  the  simplification  is  likely  to  affect,  to  a  slight  degree, 
the  calculated  pressure  distribution  mean  value  as  well  as  its  harmonic  variations. 

It  is  likely  that  the  unsteady  case  was  incorrectly  setup  for  the  PUF-10  calculations.  As  such,  this  statement  is 
not  meant  to  judge  the  capability  of  the  PUF-10  panel  method. 
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Comparisons  at  r/R  —  0.5  and  r/R  —  0.9  radii  were  also  made  for  the  mean,  first 
harmonic  amplitude  and  first  harmonic  phase.  The  comparisons  show  similar  agree- 
ment between  PUF-14  and  the  experiment.  The  comparisons  with  the  experiment 
confirm  that  PUF-14  is  correctly  solving  the  boundary  value  problem  for  this  case  of 
inclined  effective  inflow  using  modestly  skewed  propeller  blades. 

5.5.3      Propeller  AY  and  Kq 

Finally,  in  the  effort  to  compare  with  this  experiment,  many  methods  were  used 
to  validate  the  various  aspects  of  the  calculations.  Thus,  as  a  by-product  of  this 
experimental  comparison,  the  propeller  forces  and  moments  were  calculated  by  several 
different  methods.  Figure  5-19  presents  the  calculations  for  the  mean  values  of  AY 
and  Kq  at  the  design  conditions  of  the  propeller.  The  comparisons  of  the  mean  value 
of  Aj  and  Kq  suggest  that  PUF-14  is  correctly  solving  the  global  blade-row  forces 
and  moments,  but  may  be  over  predicting  the  mean  values  by  a  few  percent. 
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Figure  5-10:  Comparison  of  the  PUF-14  and  PUF-2.1  force  analysis 
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Figure  5-11:  Comparison  of  the  PUF-14  and  PUF-2.1  moment  analysis 
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Figure  5-12:  Propeller  4679  blade  shape  for  the  unsteady  calculations. 
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Figure  5-13:  Mean  pressure  distribution  on  one  blade  of  DTMB  4679  at  r/R  =  0.7  for  the  design  Js 

of  1.078. 
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Figure  5-14:  Mean  pressure  distribution  on  one  blade  of  DTMB  4679  at  r/R  =  0.7  for  Js  of  0.719. 
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Figure  5-15:  First  harmonic  amplitude  pressure  distribution  on  one  blade  of  DTMB  4679  at  r/R 
0.7  for  the  design  Js- 
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Figure  5-16:  First  harmonic  amplitude  pressure  distribution  on  one  blade  of  DTMB  4679  at  r/R 
0.7  for  Js  of  0.719. 
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Figure  5- 17:  First  harmonic  phase  pressure  distribution  on  one  blade  of  DTMB  4079  at  r/R  —  0.7 
for  the  design  J$. 
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Figure  5-18:  First  harmonic  phase  pressure  distribution  on  one  blade  of  DTMB  4079  at  r/R  =  0.7 
for  Js  of  0.719. 
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Figure  5-19:  Propeller  4679  mean  Kt  and  Kq  for  the  design  Js  =  1.078  and  an  off-design  Js  =  0.719. 
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Chapter  6 

Coupling  with  Three-Dimensional, 
Steady  Viscous  Flow  Solver 

6.1      Overview  of  Coupling 

The  concept  of  coupling  a  potential  flow  propeller  unsteady  forces  (PUF)  analysis 
method  with  a  viscous  flow  solver  is  relatively  simple.  The  difficulty  of  designing 
coupled  methods  comes  in  resolving  the  details  of  the  processes.  In  practice,  the 
coupling  of  lifting-surface  theory  with  viscous  flow  solvers  has  proven  useful  in  a.\- 
isymmetric  cases.  Some  early  efforts  of  coupling  potential  flow  analysis  methods 
with  viscous  solvers  were  by  Science  Applications  International  Corporation  [43]  and 
Kerwin,  et.al.  [28]  in  1993  and  1994.  respectively. 

Prior  to  the  implementation  of  coupling  methodologies,  perhaps  the  most  trou- 
blesome assumption  had  been  that  the  propeller  operates  in  potential  flow.  In  fact, 
propellers  almost  always  operate  in  the  highly  rotational  shear  flows  of  the  ship's 
boundary  layer  and  wake.  The  propeller  does  not  operate  in  the  flow  field  measured 
in  nominal  wake  surveys  and  one  designed  to  that  flow  would  often  be  deficient  in 
actual  use.  This  circumstance  leads  to  the  so-called  effective  wake,  or  effective  in- 
flow, problem.  The  need  for  an  effective  inflow  description  arises  from  the  flow  at 
the  propeller's  location  being  different  with  and  without  the  propeller  present.  This 
difference  is  due  to  coupling  of  the  propeller's  induced  velocity  field  to  the  vorticity 
in  the  incoming  flow.    A  redistribution  of  the  vorticity  in  the  inflow,  and  hence  the 
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shear  profile  of  the  flow  at  the  propeller,  results. 

Potential  flow  theory  provides  a  powerful  basis  for  representing  the  propeller's 
own  flow  field  but  contains  no  mechanism  to  treat  the  effective  inflow  problem.  On 
the  other  hand,  viscous  flow  methods  naturally  capture  the  vortical  phenomena  of 
the  inflow  to  the  propeller  but  offer  a  poor  framework  for  representing  the  propeller 
itself.  The  coupled  method  described  herein  couples  the  two  formulations,  viscous 
flow  and  potential  flow,  to  rationally  address  the  effective  inflow  problem  in  a  fully 
three-dimensional  flow  field. 

To  handle  the  effective  inflow  problem  requires  a  flow  solver  that  can  explicitly 
model  the  transport  of  vorticity.  There  are  a  number  of  methods  that  do  this.  An 
additional  requirement,  though,  is  that  the  solver  be  able  to  capture  separation. 
Integrated-stern  concepts  involve  hull  forms  that  actually  depend  on  the  action  of 
the  propulsor  to  inhibit  separation  as  discussed  in  Dai  et.nl.  [9].  Chen  ct.al.  [5], 
Warren  [41].  Thus,  a  viscous  flow  solver,  specifically,  a  Reynolds- Averaged  Navier- 
Stokes  (RANS)  solver  is  used. 

The  method  of  using  a  lifting-surface  with  a  viscous  solver  has  the  significant 
advantage  of  easily  supporting  multiple  blade-row  analysis.  To  treat  such  a  problem  in 
potential  flow  alone  leads  to  numerical  difficulties  as  wake-sheet  singularities  approach 
control  points  on  the  downstream  blade  row.  In  the  coupled  method,  all  the  vorticity 
is  dealt  with  by  the  RANS  code  so  that  there  are  no  singular  structures  and  the 
velocity  field  is  smooth.  In  addition,  the  potential-flow  propeller  analysis  only  has  to 
deal  with  one  blade-row  at  a  time. 

The  first  part  of  the  problem  treats  the  viscous  flow  around  a  body  while  the  second 
part  treats  the  inviscid  problem  of  the  flow  around  the  blade-row.  The  presence 
of  the  blade-row  in  the  viscous  solution  is  represented  by  a  suitable  distribution 
of  body  forces  obtained  from  the  blade-row  analysis.  In  return,  the  viscous  flow 
solver  provides  a  distribution  of  total  velocity  which  serves  as  the  input  to  the  blade- 
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row  analysis.  However,  the  precise  relationship  between  body  forces  and  blade-row 
inflow  velocity  is  not  obvious.  In  particular,  this  portion  of  the  thesis  explores  the 
extraction  of  the  effective  inflow  from  the  total  inflow  in  fully  three-dimensional  flow 
and  the  relationship  of  the  forces  to  the  total  and  effective  inflow.  Additionally,  this 
thesis  explores  the  relationship  of  the  flows  during  three-dimensional  coupling  when 
including  hub/duct  images,  multiple  blade-rows  and  a  highly  contracting  flow  field. 

Recently,  coupling  a  PUF  code  with  a  RANS  code  was  done  at  the  Naval  Surface 
Warfare  Center,  Carderock  Division  as  a  proof  of  concept  for  the  current  work  [19]. 
Their  coupling  successfully  showed  the  concept  of  three-dimensional  coupling  for  an 
open-water  single  blade-row  case.  Although,  their  coupling  used  a  three-dimensional 
RANS  code,  it  only  extracted  a  two-dimensional  inflow  for  computing  the  effective 
wake.1 

6.2     Relationship  of  Time- Average  and  Unsteady 

Imagine  a  fully-appended  underwater  body  being  towed  forward  through  the  water. 
From  the  vantage  point  of  the  body,  the  gross  flow  field  does  not  change  with  the 
advancing  of  time.  There  will  be  some  turbulent  eddies  which  seemly  fluctuate  ran- 
domly with  time,  but  the  magnitudes  of  the  fluctuations  are  small  compared  to  the 
gross  flow  field  velocity.  The  flow  field  can  be  averaged  in  time  to  remove  the  turbu- 
lent fluctuations.  Then,  from  the  vantage  point  of  the  body,  the  flow  field  is  constant 
in  time  -  or  in  the  language  of  this  thesis  -  the  flow  field  is  time-averaged. 

The  appendages  leave  a  momentum  defect  associated  with  their  wakes.  Towards 
the  body  stern,  there  are  regions  near  the  freest  ream  velocity,  regions  of  higher  velocity 
and  regions  of  lower  velocity.  By  using  process  of  time  averaging,  the  flow  field  retains 
the  flow  variations  around  the  bodv's  features. 


'The  NSWC  coupling  used  DTNS3D  and  PUF-2.  PUF-2  assumes  that  the  velocity  just  upstream  of  the  propeller 
can  be  used  everywhere  over  the  axial  extent  of  the  propeller.  Its  formulation  ignores  axial  variations  in  the  inflow 
velocity.  Thus,  the  NSWC  coupling  used  a  two-dimensional  effective  inflow  disk  at  a  predetermined  axial  location 
just  upstream  of  the  propeller  and  ignored  axial  variation  and  radial  contraction  of  the  inflow. 
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Now  imagine  the  underwater  body  has  a  propeller  attached  to  it's  stern  which 
pushes  the  body  through  the  water.  At  a  fixed  point  in  space  near  the  propeller,  the 
passage  of  the  propeller  blades  produces  a  time- varying  velocity.  As  shown  some  time 
ago  by  Breslin  [4],  this  fluctuating  velocity  field  is  due  to  a  combination  of  the  motion 
of  the  blades  relative  to  the  observer,  and  the  time-varying  loading  on  the  blades 
themselves.  This  gives  rise  to  velocity  fluctuations  at  all  harmonic  orders  of  shaft 
rate  -  not  just  blade  rate.  In  this  way,  rotating  a  propeller  through  a  time-average 
flow  field  can  produce  propeller  loading  which  varies  in  time  -  or  in  the  language  of 
this  thesis  -  produces  unsteady  propeller  loading. 

6.3     Lifting-Surface/RANS  Coupling  Details 

The  lifting-surface  method  solves  the  time-varying  boundary  value  problem.  These 
solutions  give  blade  loading  as  a  function  of  time  and  corresponding!}'  as  a  function 
of  blade  angular  position.  The  blade  loading  is  translated  into  stationary  body  forces 
in  the  RANS  domain.  The  body  force  influences  the  RANS  solution,  which  in  turn, 
influences  the  flow  field  near  of  the  propeller.  Then,  the  flow  field  is  extracted  from 
RANS  and  used  as  input  to  the  lifting-surface  solution.  The  whole  cycle  is  repeated 
until  convergence.  Figure  6-1  shows  the  coupled  solution  methodology.  Figure  6-2 
shows  the  some  of  the  details  of  an  iteration.  The  coupling  suite  of  utility  programs 
are  labeled  in  these  figures.  Their  names  are:  cs-setup  -  setup  program,  cs_unns  - 
velocity  extraction  program,  and  cs_chipr  -  force  interpolation  program. 

6.3.1      Time-Average  Induced  Velocity 

Some  means  must  be  formulated  to  get  the  effective  inflow  from  the  flow  solver.  The 
goal  is  to  extract  the  time-average  effective  inflow  from  the  time-average  total  inflow 
as  determined  by  the  viscous  flow  solver.  The  relationship  of  effective  inflow  to  the 
total  inflow  is  given  in  equation  6.1.    The  axisymmetric  effective  inflow  is  discussed 
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Figure  6-1:    Coupled  solution  methodology  flowchart 
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Figure  6-2:    Coupled  solution  iteration  details 
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in  numerous  papers  including  Kerwin  et.  al.  [28]. 

Veffective   —    *  total  ~    *  induced  (°-l) 

The  time- varying  loading  on  a  rotating  propeller  will  induce  a  time- varying  velocity 
at  a  fixed  point  in  space  near  the  propeller.  Some  means  must  be  made  to  average 
the  velocity  fluctuations  due  to  a  rotating  propeller.  Figure  6-3  shows  the  calculated 
induced  velocities  for  one  fixed  point  in  space  due  to  the  rotation  of  a  three-bladed 
propeller.  The  plotted  data  describes  the  time-varying  induced  velocity  produced  by 
all  blades  during  one  propeller  revolution.  The  figure  also  shows  a  constant  line  at  the 
value  of  the  time-average  velocity  resulting  from  averaging  the  time- varying  velocity 
produced  by  the  blade-row.  This  time-averaging  is  necessary  at  all  field  points  in  the 
swept  volume  of  the  propeller.2 

The  current  methodology  uses  the  fully  three-dimensional  swept  volume  of  the 
propeller  to  couple  with  the  viscous  solver.  As  expected,  then,  the  time-average 
induced  velocity  must  be  calculated  at  even'  control  point  at  every  blade  position. 
The  entire  swept  volume  is  described  at  field  points  which  have  been  especially  chosen 
to  coincide  with  the  control  points  whenever  the  blade  passes  the  field  point.  This 
method  removes  the  singular  effect  of  a  vortex  approaching  too  near  a  field  point. 
Additionally,  by  ensuring  the  solution  time  step,  Ng,  is  evenly  divisible  by  the  number 
of  blades,  a  field  point  is  guaranteed  to  be  co-located  with  a  control  point  at  every 
blade  time-step  position. 

After  obtaining  the  time-average  induced  velocity,  then  the  effective  inflow  can  be 
obtained  by  the  relationship  in  equation  6.1  and  seen  visually  in  figure  6-4.  Figure  6- 
4  shows  some  interesting  features.  The  upstream  flow  has  a  velocity  deficient  over 
about  SO  degrees  of  the  inflow  disk.  When  the  propeller  passes  the  region  with  the 
slower  velocity,  the  blade  loading  increases.    Higher  blade  loading  results  in  higher 

As  a  computational  note,  the  calculation  of  the  time-average  induced  velocity  is  extremely  computationally 
intensive.  The  difficulty  lies  with  the  need  to  capture  the  large  velocities  produced  by  the  blade  vortex  sheet  as  a 
field  point  is  brought  very  close  to  the  sheet.  This  need  must  be  balanced  with  the  need  to  minimize  the  required 
computational  time.  This  calculation  could  be  improved  in  the  future. 
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Velocity  at  FP  due  to  all  three  blades 


Time-Average  Velocity  at  FP 

/ 


100  200 

Blade  Position  (degrees) 


Figure  6-3:  This  specific  example  shows  the  induced  velocity  at  one  field  point  (FP)  near  mid-span 
and  mid-chord  of  the  three-bladed  propeller. 


induced  velocity  as  seen  by  the  lighter  colors  in  the  upper-right  plot.  The  velocity 
deficient  remains  in  the  effective  inflow.  Note  that  downstream  of  the  propeller,  the 
velocity  deficient  has  causes  additional  contraction  of  the  propeller  tip  streamtube 
due  the  higher  blade  loading  in  the  slower  velocity  region.  The  contraction  is  seen 
as  darker  colors  in  the  total  inflow  plot.  The  darker  colors  show  the  relatively  slower 
freestream  being  pulled  in  to  the  plotted  flow  domain. 

6.3.2      Time- Average  Forces 

It  is  desired  to  have  the  viscous  flow  solver  solve  for  the  time-average  velocity  that 
represents  the  flow  in  the  presence  of  a  propeller.  The  forces  calculated  by  the  pro- 
peller analysis  must  somehow  be  placed  into  the  flow  domain.  One  method  introduces 
the  time-average  forces  into  the  flow  domain  as  body  forces.3  In  this  way,  when  fluid 

Body  forces  are  analogous  to  gravity  forces.   Gravity  is  a  body  force  that  is  directed  downward  and  is  generally 
constant  in  strength.  A  propeller  body  force  is  placed  in  the  flow  solver  with  direction  and  strength  according  to  the 
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Total  Inflow  -  Induced  Inflow  =  Effective  inflow 


Figure  6-4:  Given  the  total  velocity  (from  RANS),  the  time-average  induced  velocity  is  subtracted 
away  leaving  the  effective  inflow  in  which  the  boundary  value  problem  is  solved. 


is  in  the  swept  volume  of  the  blade-row,  the  fluid  is  accelerated  according  to  the  forces 
produced  by  the  blade-row.  Thus,  the  fluid  feels  the  presence  of  the  blade-row  and 
the  viscous  solver  produces  a  solution  that  represents  the  time-average  total  inflow  to 
the  blade-row.  In  this  work,  based  in  part  on  the  conclusions  of  Stern  et.  al.  [37,  36], 
blade-row  forces  are  placed  in  the  viscous  flow  solver  as  body  forces  in  the  swept 
volume  of  the  blade-row. 

Time-averaging  of  the  body  forces  must  be  accomplished  such  that  the  time- 
average  body  force  is  consistent  with  the  flow-solver  time-average  velocity.  Kerwin 
et.  al.  [28]  discusses  the  relationship  of  propeller  forces  to  body  forces  for  the  circum- 


blade-row  solution. 
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ferential  mean  flow  field  as  well  as  the  relationship  of  the  total  circumferential  inflow 
to  the  effective  inflow. 

Kerwin  et.  al.  [28]  shows  that  the  forces  that  must  be  applied  in  the  flow  solver  to 
produce  the  correct  effective  inflow  arise  for  the  circumferential  mean  inflow  and  the 
blade  vorticity.  The  circumferential  mean  inflow  is  distinct  for  the  local  inflow.  The 
local  inflow  is  the  inflow  exactly  at  the  blade  mean-camber  surface.  The  circumferen- 
tial mean  inflow  is  the  inflow  that  is  averaged  at  one  field  point  during  one  propeller 
revolution. 

The  extension  to  three-dimensional  unsteady  analysis  is  as  follows.  To  place  the 
correct  forces  into  the  flow  solver,  the  forces  must  be  calculated  using  the  blade 
vorticity  and  the  time-average  velocity.  That  is.  to  keep  the  time-average  velocities  the 
same  in  the  flow  solver  and  in  the  blade  BVP,  an  equivalent  force  must  be  transmitted 
to  the  flow  solver.  The  equivalent  force  will  maintain  equal  circulations  in  the  flow 
solver  and  the  blade  BVP.  Using  a  force  calculated  with  the  blade  local  velocities 
would  be  incorrect  and  lead  to  erroneous  results  when  coupled  with  the  flow  solver. 
See  section  7.1  for  more  discussion  on  the  coupling  method,  and  reference  [28]  for 
more  discussion  on  the  forces  required  to  obtain  the  correct  axisymmetric  solution. 

6.3.3     Specific  Details  of  a  Rotor 

With  each  revolution,  the  rotor  blade  sweeps  a  constant  volume  of  fluid.  The  revo- 
lution is  discretized  by  the  number  of  time  steps  per  revolution,  No,  such  that  blade 
forces  are  calculated  at  No  locations.  Discretizing  in  this  manner  yields  the  blade 
forces  as  a  function  of  angular  position.  The  discrete  approximation  of  forces  is  made 
continuous  in  the  RANS  domain  by  interpolating  the  forces  in  the  region  between 
discrete  blade  Ng  positions.  Figure  6-5  shows  one  meridional  plane  which  contains 
rotor  forces  interpolated  onto  the  RANS  cells.  Figure  6-6  shows  the  forces  in  the 
swept  volume  of  the  rotor. 

By  interpolating  forces  between  discrete  blade  positions,  the  RANS  domain  is  pro- 
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Figure  6-5:  This  figure  shows  one  meridional-plane  containing  blade  forces  interpolated  to  the  viscous 
solver's  cells. 


vided  a  continuous  distribution  of  forces  in  the  swept  volume  of  the  rotor  outline. 
Additionally,  the  forces  in  the  swept  volume  are  time-average  forces.  Local  effects 
such  as  secondary  flows  and  tip  vortex  flows  cannot  be  captured  due  to  the  contin- 
uous nature  of  the  forces.  While  the  mean  and  unsteady  forces  are  captured,  the 
contribution  to  these  forces  due  to  secondary  flows  is  not  captured. 

6.3.4      Specific  Details  of  a  Non-rotating  Blade  Row  (Stator) 

The  stator  blades  could  be  gridded  with  RANS  cells  just  as  the  body,  duct  and 
appendages  are  gridded.  This  would  leave  only  the  rotor  to  be  modeled  with  the 
body  forces  in  the  current  methodology.  However,  the  burden  of  gridding  the  stator 
details  could  be  removed  by  modeling  the  stator  with  the  body  forces.  Given  that 
the  rotor  is  already  modeled  with  the  body  forces,  it  seems  appropriate  to  model  the 
stator  in  the  same  manner.  In  the  work  presented  herein,  the  stator  is  not  gridded 
with  the  RANS  grid. 

The  representation  of  a  stator  in  this  methodology  is  similar  to  the  rotor  rep- 
resentation with  some  unique  differences.    The  stator  is  represented  in  RANS  with 
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Figure  6-6:  This  figure  shows  a  notional  inflow  with  the  corresponding  time-average  forces  in  the 
swept  volume  of  the  rotor.  (Color  reproduction  in  figure  B-l.) 


body  forces,  which  are  only  placed  where  the  mean-camber  blade  surface  overlaps  the 
RANS  cells. 

The  correct  forces  to  represent  a  stator  in  RANS  are  those  forces  generated  with 
time-average  induced  velocities  (TAV)  just  the  same  as  for  a  rotor.  However,  for  a 
stator,  the  TAV  are,  in  fact,  identical  to  the  local  blade  induced  velocities  (LBV). 
One  unique  difference  between  the  stator  and  rotor  representation  is  that  the  stator 
forces  are  not  distributed  in  a  continuous  swept  volume  like  rotor  forces.  The  stator 
forces  are  local  to  the  blade  surface,  which  does  not  rotate.  Thus,  the  forces  are  only 
located  in  the  RANS  cells  only  where  the  stator  surface  overlaps  the  RANS  cells.  By 
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representing  the  stator  in  this  way,  the  rotor  analysis  will  capture  some  of  the  stator 
blade-rate  harmonic  forces  due  to  the  stator  interaction. 
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Chapter  7 

Preliminary  Methodology 
Validation 

7.1      Propeller  4119 
7.1.1      Setup 

The  DTMB  Propeller  4119  is  used  here  to  perform  a  preliminary  validation  of  the 
coupling  procedure.  The  case  is  close  enough  to  an  ideal  case  that  some  aspect  can 
be  verified  by  inspecting  the  results.  The  DTMB  Propeller  4119  is  a  three-bladed 
propeller  of  relatively  simple  geometry.  The  propeller  was  designed  by  Denny  for 
uniform  inflow  as  a  double  thickness  version  of  Propeller  4118  [10].  Experimental  data 
taken  by  Jessup  includes  measured  pressure  distributions,  boundary-layer  profiles, 
momentum  based  drag  coefficients  and  downstream  wake  surveys  for  this  propeller 
at  the  design  advance  coefficient  [22].  Using  the  new  coupling  methodology,  the 
measured  circulation  and  the  experimental  results  for  Kj  and  Kq  are  compared  to 
the  coupled  results 

For  the  water  tunnel  experiment,  the  diameter-based  Reynolds  number  {Reo)  was 
766,400  at  the  design  advance  coefficient  of  0.833.  The  chordwise  Reynolds  numbers 
varied  from  0.8  to  2.0  million  from  hub  to  tip.  Experimental  results  were  obtained 
with  the  blade  both  tripped  and  smooth.  The  tripping  was  done  with  roughness  at 
the  leading  edge. 
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7.1.2     Calculations 

The  propeller  is  placed  in  a  cylindrical  flow  domain  within  the  RANS  domain.  The 
chosen  RANS  grid  is  89  axially  by  81  radially  by  25  azimuthally.  The  grid  is  rather 
coarse  but  has  been  proven  adequate  for  these  results.  Figure  7-1  shows  one  merid- 
ional plane  of  the  grid  and  the  outline  formed  by  revolving  the  grid  about  the  x-axis. 
The  swept  volume  of  propeller  4119  is  shown.  The  grid  extends  4  propeller  radii 
upstream,  6  radii  downstream  and  12  radii  radially.  The  specified  upstream  inflow  is 
axisymmetric  with  U  =  1  and  V  =  W  —  0. 

Y 


Figure  7-1:  Straight-shaft  RANS  grid,  grid  outline,  and  the  swept  volume  of  propeller  4119 

It  should  be  noted  that  the  convergence  of  the  numerics  is  first  checked  using  a 
stand-alone  PUF  evaluation  in  an  assumed  effective  inflow.  The  blade  lattice  con- 
verges using  a  15x15  lattice  to  within  one  percent  of  the  converged  value  of  Kj.  The 
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K't  converges  when  using  variable  N$  within  about  a  very  coarse  10  steps  per  revo- 
lution. However.  Ng  of  10  is  unacceptable  when  performing  coupled  analysis.  In  a 
coupled. evaluation  using  PUF  and  a  RANS  code,  the  time-average  velocities  play  an 
integral  part  in  ensuring  convergence  to  the  correct  solution.  At  present,  the  calcula- 
tion of  the  time-average  velocities  uses  N$  angular  positions.  As  Ne,  becomes  smaller 
the  accuracy  of  the  calculation  for  the  time-average  velocity  becomes  much  worse. 
Therefore,  the  conclusion  is  that  Ne  should  be  a  minimum  of  about  60. 

Figure  7-2  shows  the  numerical  dynamics  during  the  coupled  analysis.  The  figure 
shows  three  different  sets  of  coupling  solutions.  The  set  marked  with  an  open  circle 
is  PBD-14  coupled  with  the  axisymmetric  version  of  DTNS.  The  set  marked  with  a 
U  is  PUF-14  coupled  with  the  three-dimensional  serial-UNCLE.  The  set  marked  with 
D  is  PUF-14  tricked  into  coupling  with  the  axisymmetric  version  of  DTNS.  Note 
this  is  an  axisymmetric  test  case.  By  using  three-dimensional  analysis  tools,  their 
preservation  of  the  axisymmetric  properties  is  check  and  the  actual  results  can  be 
directly  compared  with  the  PBD/DTNS  coupling  solution  that  has  been  verified  by 
others.  There  appears  to  be  some  error  associated  with  the  three-dimensional  coupling 
that  can  be  either  attributed  to  the  coupling  suite  of  utility  codes  or  attributed  to  the 
three-dimensional  RANS  code.  For  this  special  case,  the  three-dimensional  coupling 
should  closely  match  the  PBD/DTNS  coupling.  The  cause  of  this  three-dimensional 
error  is  not  obvious. 

During  a  coupled  analysis,  the  numerical  dynamics  followed  the  expected  trends. 
At  iteration  0  (PUF  analysis  in  the  nominal  inflow  of  U  —  1  and  V  =  W  —  0),  the  Kt 
is  expected  to  be  too  low.  The  reason  is  that  the  nominal  inflow  does  not  include  the 
induced  velocities  in  the  propeller  wake.  This  causes  the  wake  pitch  to  be  too  tightly 
coiled,  which  in  turn,  causes  too  much  induced  velocity  at  the  propeller  region.  Large 
induced  velocity  unloads  the  propeller  and  causes  Kj  to  be  too  low. 

At  the  next  iteration  (the  second  PUF  analysis),  the  Kt  is  too  high.  The  reason 
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Figure  7-2:  Coupled  analysis  convergence  history  for  propeller  4119 

is  as  follows.  Iteration  0  body  forces  are  low  which  causes  the  total  velocity  from 
RANS  to  be  smaller  that  it  should  be.  Additionally,  Iteration  0  induced  velocity, 
which  is  subtracted  from  the  total,  is  too  large  due  to  the  tightly  coiled  wake.  Both 
the  total  too  low  and  the  prior  induced  too  high  cause  the  current  effective  inflow  to 
be  too  low.  This  combination  results  in  Kj  too  large.  The  dynamics  will  eventually 
converge  to  the  correct  solution,  with  the  effective  inflow  equal  to  the  nominal  inflow 
and  with  the  propeller  wake  aligned  with  the  total  flow  field. 

One  of  the  most  useful  aspects  of  this  test  case  is  in  comparing  the  effective  inflow 
to  the  expected  effective  inflow.  The  flow  has  an  upstream  velocity  of  U  =  1  with 
V  —  W  =  0,  i.e.  circumferential-mean  inflow.  No  vorticity  exist  in  the  flow  at  the 
upstream  boundary.  The  grid  was  not  intended  to  be  fine  enough  to  accurately  cap- 
ture the  boundary  layer.  However,  the  boundary  layer,  which  does  build  on  the  shaft, 
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is  confined  to  a  small  percentage  of  the  blade  span.  Additionally,  the  experimental 
setup  uses  a  shaft  with  a  nose  cone  upstream  of  the  propeller  which  changed  the 
propeller  inflow  velocity  profiles  from  the  profile  that  is  being  calculated  here  with 
RANS.  This  effect  is  neglected  in  this  study  since  reasonable  results  can  be  obtained 
with  less  computational  time  using  a  straight-shaft  grid. 

For  this  special  case  (of  using  a  straight  shaft)  away  from  the  shaft  boundary  layer 
effects,  the  nominal  flow1  must  be  equal  to  the  effective  inflow.  This  is  true  since 
there  is  no  pre-existing  vorticity  in  the  inflow  to  interact  with  the  blade  induction. 

A  confirmation  of  the  three-dimensional  error  can  be  seen  in  the  effective  inflow. 
After  convergence  of  the  coupled  lifting-surface/RANS  solution,  the  effective  inflow 
to  the  blade  boundary  value  problem  must  be  U  —  1  with  V  =  W  =  0.  Figure  7-3 
shows  the  effective  inflow  after  the  convergence  of  the  lifting-surface/RANS  solution. 
Propeller  drag  and  thickness  have  been  set  to  zero.  This  check  must  hold  true  for  any 
propeller,  with  any  number  of  blades  and  at  any  advance  coefficient.  Many  variations 
of  the  check  were  tried  and  each  check  resulted  in  similar  effective  inflow. 

7.1.3      Possible  Source  of  Error 

For  completeness,  this  section  is  included.  Many  avenues  were  tried  to  uncover  why 
the  effective  inflow  does  not  match  the  expected  inflow.  Sources  of  error  include  a 
myriad  of  possibilities.  However,  one  source  has  stood  out.  The  two  RANS  codes, 
DTNSax  verses  UNCLE,  provide  different  velocities  when  given  the  identical  body 
force  distribution.  Note  that  the  effective  inflow  was  low  by  1-2%  when  PUF  was 
coupled  with  either  the  DTNS3D  or  UNCLE.  Leading  one  to  conclude  the  if  the  error 
is,  in  fact,  in  the  RANS  code,  then  the  error  is  in  the  body  force  routines  (which 
are  common)  or  in  a  fundamental  assumption  of  the  RANS  formulation  or  in  some 
combination.2 

That  is,  the  flow  that  would  be  present  in  the  absence  of  the  propeller. 

This  section  is  not  drawing  the  conclusion  that  the  error  is  in  the  RANS  solvers,  but  simply  stating  that  the 
difference  in  the  3-D  solvers  and  the  axisymmetric  solver  is  possibly  the  source  of  error.  More  detailed  study  of  this 
difference  in  velocity  is  necessary  before  more  definitive  conclusions  can  be  drawn. 
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Figure  7-3:  Away  from  hub  and  tip  effects,  the  axial  effective  inflow  for  propeller  4119  on  a  straight 
shaft  should  be  equal  to  1.0. 


The  3-D  RANS  solution  estimates  the  axial  velocity  to  be  a  few  percent  lower 
than  the  axisymmetric  RANS.  With  a  total  velocity  a  few  percent  lower,  the  effective 
velocity  will  be  lower.  Of  course,  the  non-linear  nature  of  coupling  a  propeller  with 
a  RANS  code  could  mitigate  or  enhance  this  effect.  But,  the  trend  it  consistent  with 
the  straight-shaft  effective  inflow  being  too  low  and  consistent  with  the  resulting 
I\'t  being  too  high.  Figure  7-4  shows  the  RANS  results  for  the  identical  body  force 
distribution.  The  body  forces  used  for  this  comparison  used  only  axial  body  forces  so 
that  any  doubts  about  coordinate  transformations  between  axi-  and  3-D  coordinate 
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systems  are  removed. 
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Figure  7-4:  Velocity  comparison  at  the  propeller  plane  (x/R  =  0)  between  DTNSax  and  UNCLE  for 
identical  axial-only  body  force  distributions. 


7.1.4      Circulation  of  Propeller  4119 

Using  Laser  Doppler  Velocimetry  (LDV)  the  flow  field  was  measured  at  two  planes 
downstream  of  the  propeller.  The  circulation  in  the  flow  field  was  computed  by 
Jessup  by  integrating  the  tangential  velocity  at  the  upstream  plane.  The  experimental 
results  are  compared  to  the  lifting-surface  calculations  in  figure  7-5.  Note  that  the 
experimental  data  shown  has  had  the  circulation  due  to  the  tangential  component  of 
the  boundary-layer  wake  subtracted.  This  was  done  by  Jessup  to  allow  comparison 
with  potential  flow  circulations. 

The  experimental  results  are  not  expected  to  exactly  match  the  calculated  results. 
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Foremost,  the  RANS  domain  does  not  include  the  nose  cone  leading  the  shaft.  Sec- 
ondly, the  differences  in  the  circulation  near  the  tip  are  due  to  the  wake  contracting 
slightly  and  a  tip  vortex  beginning  to  form  at  the  measurement  plane.  Near  the  hub, 
the  boundary-layer  calculations  may  be  under-predicting  the  flow  separation  that 
occurs  there. 

Figure  7-5  does  confirm  the  that  the  circulations  calculated  in  the  lifting-surface 
problem  are  very  near  to  the  circulations  placed  in  the  RANS  domain.  Using  theoret- 
ical extrapolations  discussed  in  section  6.3.2,  this  figure  supports  the  notion  that  the 
correct  forces  are  being  transferred  into  the  RANS  domain.  However,  this  figure  also 
indicates  that  the  coupled  lifting-surface/RANS  solution  is  resulting  in  higher  blade 
loading  than  expected,  which  can  be  reasoned  as  consistent  with  an  effective  inflow 
that  is  too  low. 

7.1.5  Blade-row  Force  Prediction 

An  additional  check  with  the  propeller  4119  is  a  check  of  the  mean  AY  and  Kq.  The 
experimental  results  at  design  J=0.833  are  AY  =  0.145  and  Kq  =  0.02S0.  Including 
the  effects  of  thickness  and  drag,  the  calculated  results  are  AY  =  0.160  and  Kq  — 
0.0293.  Recall  the  RANS  grid  does  not  accurately  model  the  nose  cone  leading 
the  shaft.  Some  error  associated  with  the  mean  AY  and  Kq  is  expected  due  to 
the  assumptions  of  the  flow  domain  geometry  and  some  error  is  associated  with  the 
incorrect  effective  inflow. 

7.1.6  Potential  Error  in  the  Coupling  Mechanics 

The  coupling  methodology  is  illustrated  in  figures  6-1  and  6-2.  The  weights  block 
shows  the  pre-calculation  for  the  interpolation.  The  weights  dictate  the  fraction  of 
forces  that  go  into  specific  RANS  cells.  Pre-calculating  weights  saves  computational 
time.  However,  if  the  weights  are  not  updated  and  the  nominal  weights  are  used 
throughout  the  coupled  analysis,  then  the  effective  inflow  will  be  in  error  by  about 


77 


Propeller  41 19 


3.5 

3 

2.5 

O      2 

o 

L1.5 

1 

0.5 

0 

-0.5 


EXPERIMENTAL  DATA  -  Vise.  Cor. 
At  x/R=0.333  in  UNCLE  solution 
PUF-1 4  Solution 


■'  '  '  I   I  I   I  I  I  '   '   I  I  '   I   I   I  I  I   '   I  '  I  I   I   I   I  *  I   I   I  I  *  I   I  I   I  I  I   I   I   '  I  I   I   u 

0.2       0.3       0.4       0.5       0.6       0.7       0.8       0.9        1 

r/R 


Figure  7-5:  Circulation  distribution  for  Propeller  4119  from  experiment,  PUF-14  lifting-surface  so- 
lution and  from  the  RANS  solution. 


2%  for  this  straight-shaft  case.  This  2%  error  in  effective  inflow  would  be  added  to 
the  error  previously  discussed.  For  this  case,  the  2%  error  in  effective  inflow  will  lead 
to  about  a  3-4%  error  in  Kj  and  Kq.  To  remove  this  possible  source  of  error,  the 
weights  are  updated  periodically  during  a  coupled  analysis. 

7.2     MIT  Pre-swirl 

7.2.1      MIT  Pre-swirl  Experimental  Setup 

In  an  attempt  to  validate  the  representation  of  the  stator  when  using  the  coupled 
methodology,  the  MIT  Pre-swirl  experiment  was  used  as  a  validation  case.  Forces 
on  the  stator  and  rotor  combination  were  obtained  experimentally  in  the  MIT  water 
tunnel.   The  stator  was  designed  by  Bowling  [3]  to  provide  swirled  inflow  to  David 
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Taylor  Model  Basin  model  propeller  4497.  The  stator  geometry  and  rotor  geometry 
are  described  in  [3].  The  experimental  geometry,  as  represented  in  the  coupled  analysis 
method,  is  shown  in  figure  7-6.  The  experiment  measured  mean  rotor  thrust  and 
torque. 


Figure  7-6:  Coupled  analysis  method  representation  of  MIT  Pre-swirl  geometry.    Flow  would  be 
generally  from  right  to  left  in  this  figure. 


The  water  tunnel  test  section  has  a  square  cross-section.  However,  as  a  simpli- 
fication, the  grid  was  made  axisymmetric  by  revolving  the  grid  around  the  X-axis. 
The  grid  outer  diameter  was  selected  so  that  the  tunnels  walls  were  represented  with 
a  cylinder  of  the  same  cross-sectional  area  as  that  of  the  test  section.  This  match- 
ing of  cross-sectional  areas  captures  the  tunnel  walls  interference  effect  as  described 
in  [13].  Figure  7-7  shows  one  meridional  plane  of  the  RANS  grid.  To  form  the  three- 
dimensional  grid,  the  meridional  plane  is  rotated  about  the  X-axis.  The  chosen  grid 
is  of  size  89  axially  by  33  radially  by  144  azimuthally.  Due  to  computational  limita- 
tions, the  grid  is  not  refined  enough  at  the  inner  and  outer  boundaries  to  capture  the 


boundary  layer. 
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7.2.2      Rotor  Experimental  Comparisons 

The  computed  force  coefficients  are  shown  in  figure  7-8.  The  forces  are  shown  ver- 
sus number  of  iterations  between  the  lifting-surface  computation  and  the  viscous 
flow  solver.  The  thrust  and  torque  computed  in  the  coupled  analysis  using  the 
PUF14/DTNS3D  method  are  within  7.3%  and  5.8%  respectively,  of  those  obtained 
experimentally.  The  same  case,  save  with  a  finer  RANS  grid,  has  been  checked  using 
PBD/DTNSax  with  thrust  and  torque  computed  within  3%)  and  1.5%  respectively,  of 
those  obtained  experimentally.  It  is  suspected  that  the  much  of  the  difference  between 
PUF/DTNS3D  and  PBD/DTNSax  results  is  due  to  the  same  effective  inflow  error 
discussed  in  section  7.1.  No  data  on  the  stator  is  available  when  operating  upstream 
of  the  rotor. 

Adjustment  to  the  calculated  rotor  thrust  due  to  the  presence  of  a  hub  vortex  was 
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Figure  7-7:  One  meridional  plane  of  the  RANS  grid  is  shown  with  the  stator  and  rotor  blades 
overlayed.  The  three-dimensional  RANS  grid  is  formed  by  revolving  this  meridional  plane  about  the 
X-axis. 
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Validation  of  PUF-14/DTNS3D  using  MIT  Pre-swirl  case 
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Figure  7-8:  Computed  MIT  Pre-swirl  force  coefficients  versus  coupling  iterations. 

calculated  by  the  method  of  Wang  [40]  as  shown  in  [7].  The  estimated  hub  vortex  axial 
force  was  less  than  0.2%  of  the  calculated  rotor  thrust  and  was  considered  negligible. 
The  change  to  the  calculated  rotor  thrust  due  to  hub  gap  forces,  as  described  in  [33], 
was  likewise  negligible. 

7.2.3     Stator  Secondary  Flows 

The  stator  is  represented  in  RANS  as  body  forces  at  the  RANS  cells  that  correspond 
to  the  stator  blade  positions  in  the  lifting-surface  code.  The  flow  field  responds  to  the 
RANS-stator  body  forces  in  a  similar  manner  as  if  a  physical  blade  was  present  in  the 
flow.  Figures  7-9  and  7-10  show  the  flow  field  at  X/ Rrotor  =  —0.62,  about  mid-chord 
of  the  stator  blade.  The  flow  accelerates  on  the  suction  side  of  the  blade  and  slows 
on  the  pressure  side.  Near  the  tip,  the  flow  tends  to  roll  over  from  the  pressure  side 
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to  the  suction  side  subsequently  forming  a  tip  vortex. 

The  major  goal  of  this  thesis  is  to  capture  the  one-per-revolution  forces  in  order  to 
predict  maneuvering  forces.  The  stator  secondary-flow  detail  will  contribute  to  the 
rotor  blade-rate  forces,  not  the  once-per-revolution  forces.  While  the  stator  secondary 
flows  are  of  interest,  they  are  not  the  center  of  effort  for  this  thesis.  Further  studies 
of  the  stator  flow  details  are  left  to  future  work. 

7.2.4      RANS-Stator  Thickness  and  Non-dimensional  Time 

The  non-dimensional  time  variable  Ng  controls  the  number  of  stator  force  positions 
written  for  conversion  into  RANS  body  forces.  The  stator  force  position  that  corre- 
spond to  stator  blade  positions  will  have  body  forces  present;  the  other  stator  force 
positions  have  forces  which  have  been  set  to  zero.3  Figure  7-11  illustrates  the  stator 
force  placement  in  RANS. 

The  RANS-stator  '"thickness"  is  defined  as  the  azimuthal  extent  of  those  RANS 
cells  where  the  forces  from  a  blade  overlap  the  RANS  cells.  Overlapped  RANS  cells 
approximate  the  stator  shape.  The  approximation  improves  with  smaller  RANS  cells 
and  higher  Ng.  The  thickness  of  the  RANS-stator  approximation  of  the  stator  blade 
will  strongly  affect  the  local  flows  around  and  behind  the  stator.  A  RANS-stator 
approximation  with  smaller  RANS  cells  leads  to  more  defined  secondary  flows  sur- 
rounding the  stator. 

When  NstatorBiad.es  =  Ng,  each  stator  must  occupy  the  entire  volume  associated 
with  its  sector  so  that  the  stator  forces  are  distributed  in  the  full  circumferential 
volume.  In  a  uniform  inflow,  the  RANS  body  forces  and  the  resulting  flow  field 
would  appear  circumferential  mean.  However,  in  the  3D  coupling,  the  apparently 
circumferential-mean  flow  field  would  not  be  an  accurate  representation  of  the  physical 
problem. 


The  program,  PUF-14,  is  written  so  that  a  stator  BVP  is  calculated  just  the  same  as  a  rotor  -  except  of  coarse, 
that  the  stator  does  not  rotate.  As  little  stator-specific  logic  as  possible  was  used  in  PUF-14.  Doing  this  may  make 
the  computation  time  longer  than  necessary  for  a  stator,  but  the  payoff  is  a  less  complicated  methodology.  Since  the 
stator  is  a  special  case  of  a  rotor,  i.e.  not  rotating,  the  variable  Ne  does  take  on  a  second  meaning. 
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Recall  that  the  forces  transmitted  to  RANS  are  the  time-average  forces.  In  a  broad 
sense,  these  time-average  forces  are  generally  obtained  as  the  product  of  the  solved 
blade  circulation  and  the  RANS  calculated  time-averaged  total  velocity.  Normally, 
a  RANS  total  velocity  is  used  to  obtain  the  effective  inflow  as  the  difference  of  total 
velocity  minus  time-average  induced  velocity.  However,  for  a  stator,  the  time-average 
induced  velocity  is  equal  to  the  local  induced  velocity.  Thus,  the  RANS  body  forces, 
which  would  be  distributed  circumferentially  when  NstatorBiades  —  Ng,  are  not  consis- 
tent with  the  local  induced  velocity.  The  corresponding  lifting-surface/R  ANS  solution 
would  be  in  error. 

The  error  due  to  inconsistency  is  removed  by  increasing  Ng  so  that  the  angular 
width  of  the  stator  body  forces  is  small.  Small  angular  width  allows  the  formulation 
to  be  more  consistent  by  causing  local  effects  to  become  more  pronounced.  Based 
on  the  MIT  Pre-swirl  case  (which  has  nine  stator  blades),  a  10°  width  distribution 
of  body  forces  in  RANS  seems  reasonable.  A  10°  RANS-stator  width  corresponds 
to  Ng  —  72  4  and  is  seen  in  figure  7-12  as  the  abscissa  of  0.1.  Interpolations  of  the 
rotor  mean  harmonic  force  convergence  plot  indicates  1%  error  in  convergence  occurs 
at  about  an  Ng  =  60.  Consequently,  a  minimum  Ng  of  60  should  be  used  for  this 
stator.  Obviously,  convergence  of  higher  harmonic  forces  requires  thinner  RANS- 
stator  width,  which  requires  Ng  to  be  higher.  For  instance,  for  this  case,  a  No  of  60 
would  yield  about  a  1%  error  in  the  mean  force  and  about  a  10%  error  in  the  9th 
harmonic  force.5  Additionally,  for  cases  with  higher  numbers  of  stator  blades,  it  is 
expected  that  thinner  RANS-stator  width  is  required  to  maintain  proper  local  effects. 
With  experience  in  stator  validation,  perhaps  the  proper  setting  for  the  RANS-stator 

360/72  =  5°,  but  the  linear  interpolation  algorithm  places  the  zero  forces  at  the  adjacent  angular  locations  so 
that  the  RANS-stator  width  becomes  2x5°  =  10°. 

5This  methodology  may  have  difficulty  if  Ng  is  very  large  and  the  RANS  cells  are  very  small  in  the  tangential 
direction.  The  difficulty  arises  because  the  RANS-stator  width  becomes  very  small  (the  stator  forces  occupy  an 
unrealistically  thin  volume  of  RANS  cells).  As  a  result,  the  stator  forces  are  introduced  into  too  small  a  volume  in  the 
RANS  domain,  which  leads  to  high  gradients  in  the  induced  velocities.  The  high  gradients  would  lead  to  artificially 
high  dissipation  near  and  downstream  of  the  stator.  This  high  dissipation  may  wash  much  of  the  stator  force  effect 
out  of  the  flow  field  prematurely.  In  the  case  of  a  stator  upstream  of  a  rotor,  the  wake  dents  from  the  upstream  stator 
may  get  too  thin  to  propagate  properly  to  the  downstream  rotor.  As  a  result,  the  higher-harmonic  rotor  forces  due 
to  the  stator's  effect  on  the  velocity  field  may  begin  to  be  incorrect. 
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width  can  be  based  on  the  ratio  of  RANS-stator  width  to  number  of  stator  blades,  or 
somehow  related  to  an  estimated  momentum  thickness. 

Physically,  the  stator  induces  very  local  effects  such  as  those  shown  in  figures  7- 
9  and  7-10.  Certainly,  it  is  possible  to  represent  the  global  influence  of  the  stator 
reasonably  well  using  circumferential-mean  lifting-surface/RANS  coupling.  To  switch 
to  circumferential-mean  stator  representation  is  a  fairly  straightforward  adjustment 
to  the  methodology.6  However,  the  present  methodology  represents  the  stator  using 
local  influences  to  capture  the  once-per-revolution  variations  in  stator  loading  which 
will  be  present  during  a  maneuver. 


To  represent  the  stator  as  a  circumferential-mean  blade  row,  two  areas  of  change  are  required.  First,  the  body 
forces  must  be  smeared  circumferentially.  For  instance,  by  adjusting  the  variable  Nq  to  equal  NstatorBlades-  Second, 
the  time-averaged  induced  velocity  must  be  calculated  using  circumferential-averaging. 
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Figure  7-9:  Local  axial  velocity  formed  by  the  stator  body  forces. 
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Figure  7-10:  In-plane  velocity  formed  by  the  stator  body  forces. 


85 


Number  of  Stator  Blades,  NB=9 

Number  of  Timesteps  /  Revolution,  Ne=72 


&> 


0© 


Ew 


Ss° 


3<e 


\o' 


,«* 


,os 


Stator  blade  forces 
%%  be  interpolated 

to  fi»  this 
RANS-stator  width 


-£ 


Figure  7-11:  RANS-stator  width  illustration 
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Figure  7-12:  Convergence  of  Kt  and  Kq  for  mean  and  9th  harmonic  rotor  forces  verses  RANS-stator 
width. 


Chapter  8 

Coupled  Analysis  of  an 
Open-Propeller  Slender  Body 

8.1     Huang  Body  1  Description 

This  test  case  is  one  of  a  series  of  axisymmetric  bodies  which  was  tested  by  Huang  with 
and  without  a  propeller  present  [18].  The  same  forebody  (DTMB  Model  5225)  was 
used  for  all  experiments  while  a  series  of  afterbodies  with  increasing  conicity  were 
examined.  Boundary  layer  profiles,  skin  friction,  pressure  tap  and  propeller  force 
measurements  were  obtained  and  used  by  Huang  to  perform  research  in  the  areas 
of  propeller/hull  interaction  and  turbulence  modeling  [16,  17,  38].  These  afterbodies 
have  been  used  extensively  for  validating  solutions  to  the  effective  wake  problem  using 
analytic  and  numerical  methods  [8,  36,  44,  45]. 

The  afterbody  considered  here,  Afterbody  1,  is  a  non-separating  stern  with  a  low 
tailcone  angle.  A  profile  view  of  the  entire  body  (Model  5225-1)  is  shown  in  figure 
8-1.  It  was  tested  in  the  presence  of  an  open  rotor  in  wind  tunnel  and  towing  tank 
facilities.  The  Reynolds  number  based  on  body  length  was  5.9  million  in  the  wind 
tunnel  tests  that  were  used  for  comparison  here. 

Propeller  model  4577  is  a  6.0-inch  diameter,  seven  bladed,  wake-adapted,  alu- 
minum propeller.  The  geometry  is  shown  pictorially  in  figure  8-2.  Open-water  per- 
formance data  for  this  propeller  was  obtained  from  propeller  4567A  using  a  propeller 
boat  in  a  towing  tank.  Propeller  4567 A  is  a  11.9  inch  diameter  duplicate  of  propeller 
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Figure  8-1:  Pictorial  representation  of  Huang  Body  1  (DTMB  Model  5225-1). 
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Figure  8-2:  Pictorial  representation  of  Propeller  4567A. 

8.2     Lifting-Surface/RANS  Coupling  Example 
8.2.1      Nominal  Flow 

The  nominal  flow  was  calculated  using  a  three-block  RANS  grid.  The  number  of  cells 
axially,  radially  and  azimuthally  was:  97x57x33,  73x57x33  and  73x57x33.  The  axial 
number,  radial  number  and  their  distribution  were  selected  based  on  convergence  of 
the  nominal  inflow.  The  number  of  cells  azimuthally  was  selected  to  provide  the 
largest  three-dimensional  grid  size  that  could  be  run  on  the  computers  at  MIT.  The 
results  for  all  yaw  angles  except  those  near  0°  are,  most  likely,  not  well  converged  with 
respect  to  the  RANS  grid;  however,  the  trends  around  0°  remain  valid  and  illustrate 
the  advantages  of  three-dimensional  coupling  using  PUF-14. 

Figure  8-3  shows  the  nominal  inflow  at  xj L  —  0.977.   The  serial-UNCLE  results 
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appear  to  over-estimate  the  axial  velocity  over  most  of  the  propeller  span.     The 
axisymmetric  boundary  layer  profiles  were  taken  from  Black  [1]. 
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Figure  8-3:  Nominal  inflow  on  Body  1  at  x/L  =  0.977  from  experiment,  axisymmetric  DTXS  and 
serial-UNCLE  models. 


Obtaining  the  Grid 

The  grid  is  obtained  using  inmesh  [6]  to  form  one  meridional  plane,  which  is  composed 
of  three  zones.  Then,  the  meridional  plane  is  rotated  about  the  X-axis  to  form  the  3D 
grid.  Figure  8-4  shows  the  overlap  of  the  blade  outline  and  the  meridional  plane.  As 
a  minimum,  at  least  10  RANS  cells  should  evenly  overlap  the  blade  outline  in  both 
spanwise  and  chordwise  directions.  The  figure  shown  here  illustrates  the  absolute 
minimum  RANS-cell/swept-blade  overlay. 

Good  overlap  of  the  blade  grid  by  the  RANS  grid  is  necessary  for  the  following 
reason.  The  lifting-surface  method  calculates  the  time-average  velocity  at  each  control 
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point.  The  calculation  is  done  using  potential  flow  analysis  and  is  therefore  exact  for 
the  calculated  blade  loading.  Large  gradients  in  the  time-average  velocities  could  exist 
across  the  blade  chord.  To  calculate  the  effective  inflow,  the  time-average  velocity 
is  subtracted  from  the  RANS  total  inflow.  If  the  RANS  flow  field  does  not  contain 
sufficient  cells  in  way  of  the  blade  to  capture  the  streamwise  velocity  gradients,  then 
the  results  from  the  subtraction  will  be  incorrect. 


"0.3  - 


Figure  8-4:  Overlay  of  the  blade  swept  outline  and  the  RANS  grid. 

8.3     Methodology  of  Coupling  with  RANS 

The  lifting-surface  methodology  uses  the  total  inflow  extracted  from  RANS  to  calcu- 
late the  effective  inflow,  then  to  solve  the  boundary  value  problem.  Figure  8-5  show 
the  converged  total  inflow  for  Huang  body  1  at  a  30°  yaw  to  port,  or  30°  drift  angle. 
To  allow  the  entire  RANS  flow  field  to  respond  to  the  propeller's  influence,  the 
time-average  lifting-surface  forces  must  be  placed  into  the  RANS  domain.  The  time- 
average  axial  force  is  shown  in  figure  8-6  at  every  3rd  position.  The  contour  variable 
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Figure  8-5:  This  figure  shows  the  total  inflow  at  every  3      velocity  disk 


is  I\'t/{Js  *  Area). 

The  time-average  velocities  provide  the  feedback  to  the  PUF  code  so  that  the  ef- 
fective inflow  can  be  determined.  By  subtracting  the  previous  time-average  velocities 
from  the  current  total  velocities  from  RANS,  PUF-14  calculates  the  effective  inflow, 
which  it  uses  to  solve  the  boundary  value  problem  at  each  timestep.  Figure  8-7  shows 
the  time-average  axial  velocity. 

8.4     Converged  Coupled  Results 

The  results  of  the  coupled  lifting-surface/RANS  method  will  produce  the  full  flow 
field  including  the  effects  of  the  propeller  on  the  flow  field.  Additionally,  the  PUF-14 
analysis  yields  the  blade-row  specifics  such  as  blade  pressure,  local-blade  forces  and 
blade-row  forces.  Figure  8-8  shows  the  convergence  of  the  mean  AY  and  Kq  at  zero 
drift  angle.1 

'Each  iteration  took  approximately  eight  hours  on  a  DEC  Alpha  333  MHz  workstation.  Of  the  eight  hours, 
about  35  minutes  were  used  to  execute  PUF-14  and  the  remainder  was  used  for  the  execution  of  serial-UNCLE.  The 
interpolation-weight  setup  program  executed  only  every  third  coupling  iteration  and  required  about  40  minutes  to 
execute. 
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Figure  8-6:  Plot  of  the  time-average  axial  force. 
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Figure  8-7:  Plot  of  the  time-average  axial  velocity  which  is  used  to  determine  the  effective  inflow. 

Table  8.1  shows  the  experimental  results  compared  to  the  calculated  results.  The 
first  four  rows  are  repeated  from  the  work  of  Black  [1].  These  rows  show  the  Propeller 
4577  performance  in  the  experimental  nominal  inflow  and  that  computed  using  differ- 
ent turbulence  models  when  coupling  PBD-14/DTNSax.  The  PUF-14/serial-UNCLE 
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results  are  shown  in  the  last  row.  The  serial-UNCLE  code  uses  a  Baldwin-Lomax  al- 
gebraic turbulence  model  which  gives  a  nominal  inflow  that  resembles  the  DTNSax 
Baldwin-Lomax  results  at  the  hub  and  the  k-e  results  from  about  mid-span  outward. 
The  results  for  the  different  turbulence  models  are  compared  in  figure  8-3.  Since 
the  nominal  inflow  most  closely  matches  the  k-e  DTNSax  results  over  much  of  the 
blade  span,  it  seems  reasonable  that  the  PUF-14  calculation  of  Kj  and  Kq  should 
also  match  the  k-e  DTNSax  results.  In  fact,  since  the  serial-UNCLE  seems  to  over 
predict  the  axial  velocity  compared  to  the  k-e,  it  is  expected  that  the  propeller  forces 
calculated  from  that  inflow  will  be  lower  than  expected  from  the  k-e  alone.  Therefore 
the  PUF-14/serial-UNCLE  results  are  expected  to  be  similar  or  slightly  lower  than 
the  k-e  DTNSax  results  for  At-  This  is,  in  fact,  the  results  that  are  obtained,  and 
this  supports  the  need  for  an  improved  turbulence  model  in  serial-UNCLE.2 


AV 

error 

KQ 

error 

Experimental  inflow 

0.2902 

- 

0.05367 

- 

PBD-14/DTNSax  k-e 

0.2520 

-13.2% 

0.04751 

-11.5% 

PBD-14/DTNSax  Baldwin-Lomax 

0.2716 

-6.4% 

0.05081 

-5.3% 

PBD-14/DTNSax  Modified  Baldwin-Lomax 

0.2876 

-0.9% 

0.05327 

-0.7% 

PUF-14/serial-UNCLE  Baldwin-Lomax 

0.2403 

-17.2% 

0.04751 

-11.5% 

Table  8.1:  Propeller  4577  performance  at  zero  drift  angle  using  various  turbulence  models. 

To  illustrate  the  prediction  of  maneuvering  forces,  the  body  was  subjected  to 
several  drift  angles.  This  simulates  that  the  body  is  in  a  yaw,  pitch  or  in  the  steady 
rotation  during  a  maneuver.  The  mean  Kj  and  Kq  are  plotted  against  the  drift  angle 
in  figure  8-9.  Note  that  due  to  computational  and  time  limitations,  only  33  RANS 
cells  were  used  in  the  azimuthal  direction  around  the  body.  Thus,  the  RANS  grid 
may  not  be  fine  enough  to  correctly  capture  the  physical  behavior  of  the  flow  field  at 
the  larger  drift  angles. 

Figure  8-10  shows  the  Huang  body  1  at  a  30°  yaw  to  port.  Streamlines  show  the 
general  direction  of  the  flow.   The  axial  propeller  force  is  shown  in  the  contour  plot 

2Note  that  the  RANS  grid  used  in  PUF-14/serial-UNCLE  coupling  is  different  from  the  grid  used  by  Black.  Slight 
differences  could  be  due  to  grid  issues. 
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on  the  blade  surfaces.  Evidence  of  the  interaction  of  the  propeller  forces  on  the  body 
RANS  solution  can  be  seen  by  the  trajectory  of  the  streamlines  passing  near  the  body 
stern. 

Figure  8-11  shows  the  trends  in  the  harmonic  amplitudes  of  the  forces  on  a  blade. 
The  forces  are  shown  in  the  xyz-coordinate  system  that  is  rotating  with  the  blade. 
Similar  trends  exist  for  the  x,  y  and  z  moments  on  a  blade.  Figure  8-12  shows  the 
plot  of  XYZ-shaft-and-bearing  forces  versus  rotation  of  the  propeller  in  the  inertial 
reference  system  fixed  to  the  body.  Similar  trends  exist  for  the  X,  Y  and  Z  moments 
on  the  shaft. 

8.5     Maneuvering  Forces 

The  coupled  lifting-surface/RANS  method  can  be  used  to  calculate  the  forces  on  a 
body  during  a  maneuver.  The  body  forces  and  moments  are  obtained  by  integrating 
the  pressure  and  drag  forces  on  the  body  from  the  RANS  solution,  and  by  adding 
the  shaft-and-bearing  forces  from  the  PUF-14  solution.  Since  the  RANS  solution 
includes  the  presence  of  the  propulsor,  the  RANS  forces  and  moments  reflect  the 
actual  thrust  required  to  move  the  body  at  the  specified  speed  (i.e.  thrust  deduction 
is  included).  The  propeller  forces  are  obtained  from  the  time-domain  lifting-surface 
solution.  Figure  8-13  shows  the  coordinate  systems  of  the  lifting-surface  analysis 
(PUF)  and  the  RANS  solution  (RANS)  relative  to  the  body-centered  coordinate 
system  (MAN)  used  for  maneuvering  force  descriptions. 

Figures  8-14  and  8-15  show  maneuvering  data  of  lateral  forces  and  yawing  moment, 
respectively.  The  "Nominal  Experiment"  is  the  measured  data  without  the  action 
of  the  propulsor.  The  experimental  nominal  data  has  a  slight  bias  which  can  be 
seen  by  the  non-zero  value  at  zero  degrees.  This  data  can  be  calculated  in  RANS 
by  converging  the  RANS  nominal  solution  (i.e.  not  coupled  with  PUF-14).  The 
calculated  data  from  RANS  alone  is  labeled  "Nominal  RANS".  The  coupled  results, 
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"Propelled  PUF-14/RANS" ,  show  the  maneuvering  forces  including  the  influence  of 
the  propeller.  No  propelled  experimental  data  was  available  for  this  case.  Note  that 
the  RANS  grid  had  only  32  cells  azimuthally  and  will  exhibit  significant  errors  when 
the  yaw  angle  becomes  large.  Calculated  data  greater  than  about  6  degrees  becomes 
erroneous  due  to  RANS  grid  azimuthal  resolution. 

The  propeller  contributes  to  maneuvering  forces  in  two  ways.  First,  the  propeller 
shaft-and-bearings  experience  forces  due  to  the  interaction  of  the  blades  with  the 
incoming  flow.  The  shaft-and-bearing  forces  directly  place  forces  on  the  body.  Second, 
the  propeller  interacts  with  the  nearby  flow  field  to  change  the  shear  and  pressure 
forces  on  the  body,  i.e.  thrust  deduction. 

Without  a  propeller  present,  the  body  in  yaw  will  experience  no  significant  heave 
force  nor  pitch  moment.  However,  with  a  propeller  operating,  both  these  components 
are  present.  Figures  8-16  and  8-17  show  this  force  and  moment,  respectively.  While 
the  magnitudes  of  these  are  smaller  than  those  shown  in  the  preceding  figures,  these 
may  have  an  important  effect  when  acting  over  a  long  time  during  a  maneuver.  The 
secondary  force  and  moment  originate  from  the  steady  shaft  force  as  seen  in  figure  8- 
12. 
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Figure  8-8:  Convergence  of  AY  and  Kq  versus  lifting-surface/RANS  iteration  for  Huang  body  1  at 
zero  drift  angle. 
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Figure  8-9:  Converged  Kj  and  Kq  versus  drift  angle  for  Huang  body  1. 
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Figure  8-10:  Huang  body  1  coupled  lifting-surface/RANS  solution  (Color  reproduction  in  figure  B-2. 
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Figure  8-11:  Force  trends  on  a  blade  in  the  blade-coordinate  system  associated  with  varying  the  drift 
angle. 
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Figure  8-12:  Force  trends  on  the  shaft  in  the  inertial  coordinate  system  associated  with  the  varying 
drift  angles. 
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Figure  8-13:  Coordinate  system  depiction 
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Figure  8-14:  Lateral  force  on  Huang  body  1 
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Figure  8-15:  Yawing  moment  on  Huang  body  1 
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Figure  8-16:  Heave  force  when  in  yaw. 
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Figure  8-17:  Pitch  moment  when  in  yaw. 
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Chapter  9 

Performance  Prediction  for  a 
Multi-Stage,  Ducted  Propulsor 

9.1      Introduction 

The  ducted  pre-swirl  unit  initially  designed  by  the  Black  [2,  28]  was  redesigned  using 
current  design  tools  during  Engineers  degree  research  by  the  author  [41].  A  description 
of  the  design  procedure  and  exact  geometry  is  presented  by  Warren  [30.  41].  This 
geometry  is  to  be  used  as  a  notional  integrated  stern  that  can  be  analyzed  by  other 
investigators.  While  experimental  data  for  this  geometry  is  not  available,  calculations 
for  this  geometry  will  be  used  to  verify  the  method  developed  herein  for  an  integrated 
stern.  A  pictorial  representation  of  the  geometry,  named  the  Sirenian,  is  shown  in 
figure  1-1. 

The  diameter  based  Reynolds  number  for  the  stator  and  rotor  were  2.3  and  1.9 
million,  respectively.  The  RANS  grid  used  here  was  derived  from  the  one  used  by 
the  author  with  modifications  to  reduce  the  computational  time  and  memory  require- 
ments of  the  serial-UNCLE  solver.  The  current  grid  has  good  discretization  axially 
and  well  as  radially.  The  y+  is  maintained  between  1  to  4  for  the  first  cell  off  the 
body.  Azimuthally,  the  grid  is  relatively  coarse,  using  only  33  cells  to  cover  the  full 
360°  arc.  As  such,  these  coupled  results  are  not  expected  to  be  accurate  at  large 
angles  of  attack.  The  azimuthal  limit  of  33  cells  was  necessary  due  to  computational 
limitations  on  the  computers  used  for  this  study. 
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It  should  also  be  noted  that  33  cells  were  chosen  to  be  an  integer  multiple  of  the  11 
stator  blades.  This  choice  reduces  harmonic  errors,  which  may  have  been  introduced 
by  the  interpolation  of  the  stator-body  forces  onto  the  RANS  grid.  Additionally, 
to  help  compensate  for  the  coarse  azimuthal  discretization,  the  stator  B-spline  was 
rotated  by  5°  about  the  x-axis.1  This  rotation  places  each  entire  stator  blade  in  one 
RANS  azimuthal  set  of  cells.  Doing  so  will  concentrate  the  body  force  and  thereby 
strengthen  the  local  velocity  gradients.  Higher  velocity  gradients  should  minimize  the 
errors  due  to  the  different  induced  velocities  calculated  in  the  RANS  domain  and  in  the 
PUF  domain.2  Finally,  the  stator  boundary  value  problem  uses  the  circumferential- 
mean  inflow.  This  is  necessitated  by  time  constraints.  Future  validation  should  use 
the  full  3-D  inflow  to  obtain  stator  side  forces  and  moments. 

9.2     Example  Case 

The  case  is  the  first  verification  of  the  whole  methodology  developed  throughout  this 
thesis  bundled  into  one  example.  The  case  is  run  at  drift  angles  of  0°,2°,4°  and  6°. 
The  total  inflow  is  interpolated  onto  a  set  up  non-radial  conical  surfaces.  The  cone 
surfaces  provide  velocity  information  in  the  region  from  just  upstream  of  the  blade- 
row  to  a  few  blade-row  radii  downstream.  Figure  9-1  shows  the  cut-away  volume 
representation  of  the  cone  surfaces  for  the  stator  in  4°  drift  angle.  The  stator  blade 
outlines  are  shown  in  the  figure. 

The  body  stern  is  shown  in  figure  9-2  with  a  cut-away  view  of  the  duct.  The  stator 
and  rotor  transition  wakes  are  also  shown.  The  wake  shapes  conform  to  the  body 
geometry,  subject  to  boundary  layer  interactions.  The  wakes  do  not  follow  cylindrical 
geometry  assumptions  that  are  present  in  previous  unsteady  lifting-surface  methods. 
These  examples  of  the  wakes  illustrate  the  benefit  of  wake-adapted  modeling.  Addi- 
tionally, the  wake  interaction  does  not  take  place  using  potential  flow  singularities; 

'Alternately,  the  RANS  grid  could  have  been  rotated. 

2The  stator  B-spline  rotation  is   not  necessary,  but  should  partly  compensate  for  the  coarse  RANS  azimuthal 
representation.  If  the  azimuthal  number  of  cells  was  much  greater,  say  110,  then  this  rotation  would  not  be  necessary. 
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Figure  9-1:  This  figure  shows  the  volume  used  to  perform  the  stator  analysis.  This  case  is  at  4°  drift 
angle.  Tangential  velocity  is  shown  in  the  contour  (Color  reproduction  in  figure  B-3). 


instead  the  body  force  modeling  in  the  viscous  solution  accounts  for  the  interaction. 
The  stator  blade-row  is  analyzed  using  potential  flow  techniques  in  a  specified  inflow 
such  as  figure  9-1,  then,  separately,  the  rotor  is  analyzed  in  its  specified  inflow.  Thus, 
each  blade-row  analysis  is  performed  independently  and  singularities  from  different 
blade- rows  do  not  explicitly  interact.  Both  the  stator  wake  and  the  rotor  wake  follow 
the  circumferential-mean  total  flow  field.  On  average  the  wakes  should  be  force  free; 
however,  no  explicit  effort  enforces  them  to  be  force  free.  Of  course,  the  assumption 
of  a  force  free  wake  also  neglects  secondary  effects  such  as  wake  defects  from  blade 
boundary  layers,  cascading  blockage  effects  from  downstream  blade-rows,  etc.  and  is 
subject  to  the  assumption  of  time-average  viscous  modeling. 

The  serial-UNCLE  forces  and  moments  are  calculated  by  integrating  the  pressure 
and  shear  stress  on  the  body  and  duct.  The  stator  and  rotor  forces  and  moments  are 
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Stator  and  Wake 


Figure  9-2:  Wake-adapted  grids  interact  through  the  viscous  solver. 

calculated  by  the  lifting-surface  method.  Figure  9-3  shows  the  lateral  force  density 
contours  on  the  body,  duct  and  blade-rows.  The  body  contours  are  the  lateral  force 
in  the  RANS  coordinate  system  shown  in  the  figure.  The  contours  on  the  blades  are 
in  the  blade-centered  coordinate  system,  which  rotates  with  the  blade  (i.e.  contours 
are  the  tangential  force).  All  contours  are  non-dimensionalized  by  the  body  radii  to 
be  consistent  with  the  serial-UNCLE  output.  This  figure  illustrates  the  lateral  force 
density  for  a  body  at  a  4°  drift  angle  from  the  starboard  bow.  The  stagnation  point 
on  the  body  has  shifted  to  starboard  and,  in  general,  the  starboard  side  of  the  body 
feels  a  force  directed  to  port.  The  stern  shoulder  has  a  strong  force  to  starboard  that 
is  generated  as  the  flow  accelerates  at  the  shoulder  due  to  both  the  stern  contraction 
and  the  yawed  inflow. 

The  flow  that  is  a  few  propeller  radii  away  from  the  propulsor  is  not  strongly 
influence  by  the  propulsor.  However,  in  the  propulsor  region,  the  forces  on  the  body 
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are  influenced  by  the  stator  and  rotor  interactions.  The  stator  body  forces  are  at 
discrete  locations  in  the  RANS  domain,  admittedly  with  a  somewhat  broad  angular 
representation.  At  the  stator  root,  the  triangular  contour  patches  of  higher  lateral 
forces  indicate  the  discrete  stator  interactions.  That  is,  forces  are  present  where  the 
stator  blade  overlaps  the  RANS  cells  and  no  forces  are  present  in  the  void  in  between 
the  stator  blade.  The  rotor  forces,  being  time-averaged,  fill  the  entire  swept  volume 
of  the  rotor  blade-row  with  forces  that  smoothly  vary  in  all  spatial  coordinates.  The 
rotor  body  forces  influence  the  RANS  solution  with  an  axial  pressure  rise.  The  rotor 
influence  is  most  clearly  seen  by  examining  the  inside  of  the  duct.  On  the  inner  port 
side  on  the  duct,  the  lateral  force  varies  from  strong  negative  to  strong  positive.  This 
arises  due  to  the  pressure  drop  due  to  the  stator  and  then  a  pressure  rise  due  to  the 
rotor. 
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Figure  9-3:  Force  contours  on  the  body,  duct  and  both  blade-rows  (Color  reproduction  in  figure  B-4.) 
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During  the  course  of  this  work  and  other  related  work,  one  major  concern  was 
noted  associated  with  the  viscous  modeling.  For  some  cases,  such  as  this  Sirenian 
body  but  not  the  Huang  Body  1  present  in  chapter  8,  the  RANS  viscous  solution 
oscillates  around  a  mean  value.  The  amplitude  of  oscillation  does  not  decrease  with 
more  RANS  iterations  and  the  residuals  do  not  improve. 

The  converged-coupled  solution  for  Sirenian  at  4°  drift  angle  was  run  until  the  ro- 
tor thrust  and  torque  converged.  The  RANS  solution  had  also  converged  as  indicated 
by  no  further  reduction  in  the  residuals.  However,  the  RANS  solution  had  not  actu- 
ally converged  to  a  stable  solution,  instead  the  RANS  solution  was  oscillating  within 
a  narrow  range.  More  iterations  of  the  RANS  solution  did  not  stabilize  the  solution; 
more  iterations  continued  to  oscillate  the  solution  about  some  mean  value.  Presum- 
ably, the  center  of  oscillation  is  the  correct  solution.  The  dotted  lines  in  figure  9-4 
show  the  lateral  force  integrated  around  each  body  section  from  RANS  solutions  that 
were  several  hundred  iterations  apart.  The  bold  line  in  an  average  of  all  solutions 
spanning  between  the  dotted  lines.  The  dotted  lines  show  some  typical  extremes  in 
the  variations. 

In  addition  to  the  oscillation  of  the  solution,  it  is  surprising  that  the  lateral  force 
integration  does  not  vary  smoothly.  For  this  Sirenian  body,  the  geometry  varies 
smoothly  from  the  bow  to  the  stern.  The  duct  is  the  only  discontinuity.  Presumably, 
if  many  more  force  predictions  (dotted  lines)  were  averaged,  then  the  average  lateral 
force  (solid  line)  would  be  smoother.  It  is  troubling  to  find  the  solution  by  averaging 
many  "solutions'1.  Interestingly,  the  force  integration  over  the  body  stern  does  not 
seem  to  oscillate.  The  observed  oscillation  could  simply  be  due  to  poor  grid  quality, 
or  could  indicate  a  more  fundamental  problem.  Similar  oscillations  have  been  noted 
in  axisymmetric  RANS  solvers  on  other  body  shapes.  Admittedly,  in  the  hands  of  a 
more  accomplished  RANS  user,  these  oscillations  may  not  be  present.  This  concern 
is  presented  for  completeness  and  is  not  meant  to  invalidate  the  viscous  modeling  in 
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RANS,  but  to  simply  point  out  an  area  of  concern.  It  is  recognized  and  supported 
that  viscous  modeling  in  RANS  offers  an  extremely  powerful  tool  for  computational 
modeling.  However,  since  the  effort  of  this  thesis  in  not  directed  toward  the  RANS 
modeling,  research  in  the  RANS  modeling  is  left  to  other  universities  and  research 
centers.  This  thesis  effort  centers  on  the  lifting-surface  blade-row  modeling  and  the 
interaction  with  the  given  RANS  solution. 


Body  Average 


Duct  Inner  Surface 


Duct  Outer  Surface 


Figure  9-4:  Lateral  force  [Fz)  verses  body  position 

9.3     Maneuvering  Forces 

Figures  9-5  and  9-6  show  calculated  maneuvering  data  of  lateral  forces  and  yawing 
moment,  respectively.  These  forces  and  moments  are  non-dimensionalized  by  the 
body  length.  No  experimental  data  are  available.  The  trends  follow  those  observed 
in  the  Huang  Body  1  validation  in  chapter  8.  The  calculated  data  from  RANS  alone 
is  labeled  "Nominal  RANS".  The  nominal  case  is  expected  to  be  on  the  verge  of 
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separating  at  the  body  stern  due  to  the  relatively  full  stern.  While  no  evidence  of 
separation  was  observed  in  the  nominal  solution,  it  remains  important  that  the  viscous 
solver  be  able  to  capture  separation.  At  larger  yaw  angles,  separation  may  be  more 
likely.  The  calculated-coupled  results,  labeled  "Propelled  PUF-14/RANS",  show  the 
maneuvering  forces  including  the  influence  of  the  propeller.  Note  that  the  RANS  grid 
used  in  this  analysis  had  only  33  cells  azimuthally  and  may  exhibit  significant  errors 
the  at  larger  yaw  angles. 

The  propulsor  contributes  to  maneuvering  forces  in  several  ways.  First,  the  rotor 
shaft-and-bearings  experience  forces  due  to  the  interaction  of  the  blades  with  the 
incoming  flow.  The  shaft-and-bearing  forces  directly  place  forces  on  the  body.  Second, 
the  stator  experiences  forces  that  act  directly  on  the  body.  Recall  that,  due  to  time 
constraints,  this  stator  analysis  is  performed  in  the  circumferential-mean  inflow,  which 
results  in  zero  side  forces.  Still  another  contribution  from  the  propulsor  is  a  change 
to  the  nearby  flow  field  which  changes  the  shear  and  pressure  forces  on  the  body.  i.e. 
thrust  deduction.  Final,  the  duct  interacts  with  each  of  these  modify  both  the  forces 
on  the  components  and  on  itself.  The  coupled  lifting-surface/RANS  method  provides 
the  mechanism  to  quantify  their  interaction.  In  this  example,  the  propulsor  effects 
on  the  maneuvering  forces  are  relatively  small,  approximately  7-10%  change  from  the 
nominal  solution. 

Figures  9-7  and  9-8  show  the  heave  force  and  pitch  moment,  respectively.  The 
secondary  force  and  moment  originate  from  the  steady  shaft  force  and  moment.  In 
both  these  figures,  some  RANS  error  can  be  seen  by  observing  that  the  nominal 
calculated  results  should  have  zero  heave  force  and  pitch  moment.  Like  the  Huang 
Body  1  case,  the  magnitudes  of  the  heave  force  and  pitch  moment  are  smaller  than 
lateral  force  and  yaw  moment,  respectively.  While  the  heave  force  and  pitch  moment 
are  smaller,  these  may  have  an  important  effect  when  acting  over  a  long  time  during  a 
maneuver.  Additionally,  the  propulsor-induced  force  and  moment  strongly  depend  on 


111 


the  body-boundary  layer  interaction  with  the  propulsor  and  may  be  relatively  small 
for  this  example. 
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Figure  9-5:  Lateral  Force 
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Figure  9-6:  Yawing  Moment 
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Figure  9-7:  Heave  force  when  in  yaw. 
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Figure  9-8:  Pitch  moment  when  in  yaw. 
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Chapter  10 

Conclusion  and  Recommendations 


An  original  computer  code,  PUF-14,  has  been  written  to  support  the  new  methodol- 
ogy developed  throughout  this  thesis.  The  lifting-surface  method,  PUF-14,  is  coupled 
with  a  RANS  code  provided  by  our  ONR  sponsor.  The  RANS  code  is  used  to  obtain 
the  inflow  velocity  field  that  the  lifting-surface  code  uses  to  calculate  the  unsteady 
forces.  Unsteady  forces  are  generated  due  the  rotating  propeller  in  a  spatially-varying 
inflow.  Time-averaged,  but  spatially-varying  body  forces  are  introduced  into  a  three- 
dimensional  volume  to  represent  propulsor  stages  in  the  RANS  flow  field.  The  entire 
RANS  flow  field  responds  to  the  blade-row  presence.  In  turn,  the  RANS  flow  field  is 
used  again  for  the  lifting-surface  analysis  of  the  blade-row.  By  alternately  updating 
the  lifting-surface  and  the  RANS  solutions,  the  blade-row  forces  and  RANS  flow  field 
converge  to  the  appropriate  solution. 

The  new  treatment  of  unsteady  force  calculations  should  greatly  improve  propul- 
sor prediction  capabilities.  The  new  treatment  is  believed  to  be  practical  in  both 
computational  load  and  in  representing  the  physical  hydrodynamic  characteristics  of 
today's  complex  propulsors. 

10.1     Improvements  on  the  Current  Method 

In  developing  this  methodology,  many  possible  improvements  have  become  obvious 
and  should  be  implemented  given  time  and  resources.  Some  improvements  may  make 
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only  small  improvements  in  accuracy;  others  may  be  more  significant.  The  improve- 
ments are  detailed  below  in  no  particular  order. 

•  Wake  and  blade  lattice  model  which  includes  rotational  variation 

The  current  method  assumes  the  streamtubes  which  convect  the  shed  and  trailing 
vorticity  are  the  circumferential-mean  streamtubes.  This  assumption  allows  one 
set  of  influence  functions  to  be  used  at  all  blade  positions.  An  improvement  would 
model  the  rotational  variations  in  the  wake  lattice  so  that  the  wake  follows  the 
full  three-dimensional  velocities. 

•  Stand-alone  unsteady  wake  alignment 

Currently,  to  align  the  wake  when  not  coupling  with  RANS,  convection  velocities 
are  used  to  stretch  the  wake  to  account  for  induced  velocity.  Therefore,  it  may 
not  be  free  of  forces  as  it  should  be.  A  wake  alignment  routine  would  enhance 
the  use  of  PUF-14  when  used  as  a  stand-alone  analysis  tool. 

•  Trailing  edge  sheet  model 

More  research  is  needed  to  extend  the  2-D  Lagrange  Kutta  condition  to  3-D. 
The  2-D  model  seems  very  promising  and,  once  expanded  into  3-D.  may  greatly 
improve  the  accuracy  of  higher  harmonic  forces. 

•  Tip  gap  model 

The  current  methodology  does  not  attempt  to  account  for  real-fluid  effects  at 
the  blade  tips.  A  semi-porous  tip  gap  model  could  be  incorporated  into  this 
method.  The  model  could  be  based  on  the  circumferential-mean  tip-gap  effects 
or  could  encompass  a  time-varying  nature. 

•  Viscous-load  and  thickness-load  coupling 

The  growth  of  viscous  boundary  layers  on  the  propulsor  blades  result  in  changes 
to  the  inviscid  pressure  distribution  and  generally  a  loss  in  lift.  The  displacement 
thickness  of  the  boundary  layer  effectively  changes  the  camber,  thickness  and 
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angle  of  attack  of  each  blade  section.  Simularily,  the  thickness  of  the  blades 
influence  the  blade  boundary  layer.  The  effects  become  more  important  for 
advanced-blade  sections.  The  time-varying  nature  of  the  viscous  effects  would 
need  consideration. 

•  Blockage 

Blade  thickness  results  in  a  blockage  to  the  flow  which  affects  the  flow  distribution 
past  the  propulsor  and  the  duct,  if  one  is  present.  The  creation  of  the  viscous 
boundary  layer  on  the  blade  surface  creates  additional  velocities  that  will  have 
an  effect  on  the  through-flow  and  performance  of  downstream  stages.  These 
factors  need  to  be  considered  if  more  accurate  propulsor  performance  is  to  be 
predicted. 

•  Improved  algorithm  for  time-average  induced  velocities 

By  far,  the  most  calculationally  intensive  portion  of  the  lifting-surface  method  is 
the  calculation  of  the  time-average  induced  velocities.  This  calculation  requires 
the  effects  at  every  control  point  during  one  blade-interval  passage  from  every 
vortex  element  on  every  blade  and  every  sheet.  A  more  advanced  algorithm 
would  improve  the  accuracy  and  greatly  improve  the  speed  of  this  routine. 

•  Coupling  with  RANS  using  a  single  disk  of  velocities 

The  guiding  cases  used  in  developing  this  methodology  were  multi-component 
propulsors  with  highly  contracting  stern  flows.  While  this  method  works  fine 
for  open-water  propellers,  the  additional  computational  load  is  a  burden.  Us- 
ability would  be  greatly  improved  for  open-water  propellers  if  the  method  were 
reduced  to  coupling  with  RANS  via  a  single  disk  of  velocities  upstream  of  the 
propeller  plane  instead  of  the  full  3-D  volume  of  velocities.  Since  RANS  veloci- 
ties would  not  be  available  downstream  of  the  propeller  to  convect  the  wake,  a 
wake  alignment  technique  would  be  required  to  ensure  the  wake  remains  force 
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free. 

•  Improved  algorithm  for  tracking  2-D  streamlines 

The  current  algorithm  for  tracing  the  2-D  streamlines  in  the  circumferential- 
mean  flow  is  not  as  robust  as  it  could  be.  The  algorithm  struggles  in  the  viscous 
sublayers  and  leads  to  difficulties  when  automating  the  coupling  process  for 
certain  geometries. 

•  RANS  turbulence  modeling 

The  serial-UNCLE  used  for  the  latter  work  implemented  a  Baldwin-Lomax  tur- 
bulence model.  This  was  proven  to  be  inadequate  in  the  Huang  body  1  compar- 
ison. More  research  should  be  directed  to  improve  the  turbulence  models  for  the 
axisymmetric  boundary  layer  flow  and  for  the  vortex-core  flows  associated  with 
maneuvering  bodies. 

•  RANS  parallel  processing 

There  is  current  research  in  RANS  parallel  processing  algorithms.  The  method 
developed  herein  should  be  coupled  to  a  parallel  version  of  a  RANS  code.  By 
far,  the  component  in  the  coupled  lifting-surfaee/RANS  method  requiring  the 
most  computer  time  and  computer  memory  is  the  RANS  code.  Effort  should 
continue  to  be  directed  at  improving  this  component. 

10.2     Possibilities 

The  stand-alone  lifting-surface  method  may  be  useful  in  other  applications.  With 
some  modifications,  one  possible  use  could  be  to  evaluate  the  transitory  response 
to  step  inflow  changes.  Both  the  transitory  blade  forces  and  the  shaft-and-bearing 
transitory  forces  could  be  determined.  Transitory  response  may  be  useful  in  designing 
the  motor  controller  for  a  propeller. 

The  coupled  lifting-surface/RANS  method  could  be  used  to  model  the  appendages 


118 


of  a  body.  For  example,  a  control  surface  is  not  too  different  from  a  stator  blade.  Thus, 
the  appendages  could  be  modeled  as  "stator"  blades;  A  four  bladed  "stator1'  would 
model  four  stern  control  surfaces.  If  the  four  blades  were  not  oriented  identically 
about  the  body  centerline,  then  perhaps  they  could  be  modeled  as  four  one-bladed 
stators.  While  the  flow  details  would  not  be  precisely  correct,  for  a  body  at  an  angle 
of  attack,  the  gross  flow  would  be  correctly  captured.  The  coupled  method  would 
provide  a  relatively  quick  estimate  of  the  maneuvering  forces  on  an  appended  body, 
without  ever  having  to  remake  a  RANS  grid. 

Finally,  another  possibility  might  use  the  coupled  method  interactively  with  a  time- 
domain  maneuvering  simulation.  The  propulsor  forces  and  moments  could  be  calcu- 
lated at  discrete  time  steps  in  a  maneuvering  simulation.  The  resulting  body-forces 
could  update  the  RANS  solution,  while  the  propulsor  forces  influence  the  trajectory 
of  the  body  in  the  RANS  domain. 

10.3     In  Retrospect 

First,  it  must  be  acknowledged  that  neither  PUF-14,  the  coupling  methodology  nor 
RANS  codes  are  perfect.  They  all  have  shortcomings  and  all  require  careful  attention 
to  avoid  the  "garbage  in  -  garbage  out"  syndrome.  With  effort,  these  methods  can 
improve. 

The  methodologies  developed  and  incorporated  into  the  stand-alone  PUF-14  pro- 
vide the  modern  propulsor  designer  a  tool  to  analyze  trends  in  propulsor  perfor- 
mance. Without  RANS,  the  stand-alone  methodology  suffers  by  not  possessing  an 
automated  wake  alignment  method.  User  input  convection  velocities  fill  this  short- 
coming. However,  as  a  propulsor  analysis  method,  the  less  user  input  to  the  physical 
hydrodynamics  the  better  the  method. 

The  coupled  methodology  consists  of  PUF-14  and  a  RANS  solver.  This  method 
has  great  potential  and  should  be  invaluable  to  the  modern  propulsor  designer.    In 
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addition  to  trend  analysis,  the  coupled  methodology  should  provide  relatively  good 
agreement  with  experiment.  The  single-most  important  advantage  is  the  ability  to 
discretely  model  the  component  stages  while  avoiding  the  complex  and  computa- 
tionally intensive  modeling  of  the  rotating  blades  in  the  RANS  domain.  It  is  my 
sincerest  hope  that  the  methodology  developed  herein  becomes  invaluable  to  today's 
and  tomorrow's  propulsor  designers. 
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Appendix  A 

A  Numerical  Method  for  the 
Computation  of  Unsteady 
Propeller  Forces 

Introduction 

During  the  computation  of  unsteady  forces,  many  possible  approaches  could  be  adopted 
in  locating  the  shed  vortices  within  their  corresponding  time  interval.    Nearly  any 
placement  will  converge  to  the  correct  answer  given  sufficiently  small  time  steps.  Ob- 
viously, the  best  numerical  scheme  is  one  that  converges  within  the  desired  accuracy 
with  minimal  computations,  i.e.  the  largest  possible  time  step. 

In  general,  the  spacing  in  the  wake  must  be  on  the  same  order  as  the  last  spacing 
on  the  foil  to  get  reasonable  results  near  the  trailing  edge.  Fine  spacing  near  the 
trailing  edge  leads  to  the  necessity  of  similarly  fine  spacing  in  the  wake  to  remove 
the  singular  behavior  at  the  trailing  edge.  Such  fine  spacing  in  the  wake  leads  to 
considerable  computational  effort. 

Past  unsteady  vortex-lattice  computational  methods  have  used  constant  spacing 
on  the  foil  with  similar-sized  spacing  extended  into  the  wake.  Spacing  arranged  in 
this  manner  has  provided  good  accuracy  of  the  global  forces  with  reasonable  sized 
elements  in  the  wake. 

However,  it  is  desired  to  resolve  the  gradients  better  on  the  foil,  especially  near 
the  leading  and  trailing  edges.    A  natural  solution  is  to  use  Glauert  constant  angle 
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Figure  A-l:  Discrete  vortices  and  control  points  along  the  chord. 

spacing,  alias  cosine  spacing,  on  the  foil.  This  choice  places  a  fine  grid  of  points  near 
the  leading  and  trailing  edges,  and  a  relatively  coarse  grid  in  the  mid-chord  region. 
Figure  A-l  compares  constant  and  cosine  spacing  on  a  two-dimensional  foil. 

Trailing  Edge  Singularity  with  Cosine  Spacing  on  the  Foil 

As  shown  in  figure  A-2,  the  flat-plate  solution  in  a  gust  is  sensitive  to  the  placement 
of  the  vortex  with  the  time-step  element.  Numerical  studies  such  as  Frydenlund  and 
Kerwin  [12]  have  examined  the  sheet  strength  as  it  transitions  into  the  wake.  The 
sheet  strength  should  be  continuous  and  smooth,  except  for  a  slope  change  at  the 
trailing  edge. 

The  location  of  the  discrete  vortex  within  the  first  time  step  element  in  the  wake 
causes  dramatic  changes  in  the  sheet  strength  on  the  foil  and  the  wake.    A  closer 
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Figure  A-2:  Trailing  edge  behavior  of  the  vortex  sheet  strength. 

examination  indicates  that  using  an  implicit  Kutta  condition  allows  the  last  control 
point  to  be  in  close  proximity  to  the  vortex  in  the  wake.  The  closeness  was  avoided 
in  previous  formulations  by  using  constant  spacing  on  the  foil  and  an  explicit  Kutta 
condition. 

Figure  A-2  illustrates  the  strong  dependence  of  the  sheet  strengths  on  the  posi- 
tion of  the  first  wake  vortex  and  concisely  illustrates  the  need  for  this  formulation. 
The  figure  describes  a  linearized  flat-plate  in  a  sinusoidal  gust  at  an  arbitrary  time 
during  the  cycle.  The  case  with  constant  spacing  on  the  foil  serves  as  a  guide  to  the 
correct  behavior.  All  cosine-spaced  cases  have  the  identical  spacing  on  the  foil  and 
wake  except  that  the  vortex  placement  within  the  wake  time  step  element  is  varied. 
Although  far  from  correct,  placing  the  vortex  at  the  quarter  chord  location  within 
the  element  seems  to  provide  the  most  reasonable  compromise.  Perhaps  not  quite  so 
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coincidentally,  the  quarter  chord  placement  of  the  vortex  along  with  the  three-quarter 
chord  placement  of  the  control  point  yields  the  exact  solution  in  steady  flow  in  the 
limit  of  an  infinite  number  of  vortices.  [20] 

This  appendix  discusses  the  formulation  and  implementation  of  a  method  to  over- 
come the  singular  behavior  of  a  discrete  wake  vortex  near  the  foil  control  point.  The 
formulation  drastically  reduces  the  computational  effort  while  obtaining  very  similar 
accuracy. 

Formulation  of  the  Two-Dimensional  Problem 

Consider  a  two-dimensional  thin  hydrofoil  advancing  with  constant  speed  U,  which 
may  be  passing  through  a  spatially  varying  velocity  field.  The  linearized  problem  will 
be  solved.  A  flat  plat,  with  chord  length,  C,  is  situated  on,  or  close  to  the  interval 
(0,C)  of  the  x-axis  in  a  flow  with  speed  U  oriented  in  the  positive  x  direction. 

The  integral  equation  for  the  distribution  of  vorticity  *){x,t)  over  the  foil  may  be 
written  in  the  following  form. 


v(x,t)+  /     -de,  -f  /      dr]  =  0  (A.l 

Jo     x  —  c  Jc       X  —  77 


where  £  and  77  are  dummy  coordinates  on  the  foil  and  wake,  respectively,  and  -)u,  is 
the  strength  of  the  shed  vorticity  in  the  wake. 

According,  the  Kelvin's  theorem,  the  total  circulation  of  the  system  must  remain 
zero, 

/    7UVMC+  /     7wOM)«*7  =  0  (A-2) 

Jo  Jc 

so  that  the  change  in  circulation  on  the  foil  must  be  followed  by  an  equal  change  in 
circulation  in  the  wake. 

The  formulation  is  completed  with  a  statement  of  the  Kutta  condition,  which  re- 
quires that  7(2:,  t)  be  bounded  at  the  trailing  edge.  The  Kutta  condition  can  be  made 
explicit  by  fixing  the  last  bound  vortex  strength  to  be  a  value  which  satisfied  the 
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desired  behavior.  Conversely,  the  Kutta  condition  can  be  implicit  which  is  accom- 
plished by  relative  placement  of  the  last  control  point  on  or  very  near  to  the  trailing 
edge. 

Discrete  Time  Step  Solution  of  the  Boundary  Value  Problem 

For  the  discrete  vortex  model  employed  here,  the  governing  integral  equation,  eq.  A.l, 
is  reduced  to  a  system  of  linear  algebraic  equations: 

£(r*);*u  +  Br5),w/-,,  =  -v;  (A.3) 

j=i  fc=i 

The  quantities  are  as  follows: 

•  Ar  is  the  number  of  bound  vortices 

•  Nw  is  the  number  of  vortices  in  the  wake 

•  BUJ  are  the  influence  functions  which  describe  the  induced  velocity  on  control 
point  i  due  to  a  unit  strength  vortex  j  located  on  the  foil 

•  W{tk  are  the  influence  functions  which  describe  the  induced  velocity  on  control 
point  i  due  to  a  unit  strength  vortex  k  located  in  the  wake 

•  \ ~'i  are  the  velocities  acting  at  control  point  i  due  to  the  boundary  conditions 

•  (T^j-are  the  unknown  bound  vortices  at  position  j  on  the  foil  at  the  current  time 
step 

•  (Ts)k  are  the  shed  vortices  at  position  k  in  the  wake. 

The  notation,  Tu6m,  represents  the  discrete  vortex  strength  in  the  time  step  el- 
ement. The  superscript  u  in  Tusm  indicates  that  the  circulation  is  counted  in  the 
un-subdivided  intervals,  which  are  counted  with  the  index  m.  The  numerical  objec- 
tive of  this  model  is  to  represent  the  rus0  with  an  alternate  set  of  shed  vortices,  Fsk 
where  k  =  1,2,-  ••,  Nj.  The  vortices  rusm  for  m .  >  0  are  left  untouched  and  equal 
Tsk  for  k  >  Nj. 

The  total  circulation  shed  into  the  most  recent  time  step  element  m  =  0  in  the 
wake  is  designated  as  Tus0.  The  Tn  represents  the  sum  of  the  bound  vortex  strengths 
at  time  step  n.  Thus,  Kelvin's  theorem  can  be  discretized  in  equation  A. 4  by  writing: 

rus0  =  -rn  +  rn_!  (A.4) 
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Profile  View  of  a  2-D  Foil 
and  its  wake 
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Figure  A-3:  Pictorial  of  subdivisions  in  the  first  time  element 

r»  =  {S(r0iJ  (A.5) 

This  new  formulation  of  the  linearized  two-dimensional  foil  differs  from  previous  for- 
mulations in  that  the  most  recent  time  step  element  is  subdivided  into  Nj  intervals. 
Each  interval  is  represented  by  a  discrete  vortex.  Thus,  the  shed  vortices  are  sepa- 
rated into  those  within  the  most  recent  time  step,  which  are  unknown  quantities,  and 
into  those  from  prior  time  step  solutions  which  are  known.  Rewriting  equation  A. 3, 
we  get: 

T,^b)iB{J  +  ,£(Ts)kWijg  =  -Vi-     £    (Ts)kWitk  (A.6) 

i=i  fc=i  k=Nf+i 

where  (Ts)jt  are  known  for  k  >  Nj  but  unknown  for  k  <=  Nj. 
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In  order  to  maintain  finite  loading  at  the  trailing  edge  of  the  foil,  the  vortex  sheet 
must  be  continuous  from  the  foil  into  the  wake.  Thus,  by  the  relationship 

*  =  -r!<r('-F»  (A-7) 

we  would  like  for  the  discrete  vortex  to  be  at  least  of  order  two  so  the  sheet  will 
be  continuous.  One  simple  way  is  to  assume  the  T(t)  is  a  polynomial  of  order  A 
and  use  Lagrange  interpolation  to  represent  the  discrete  strengths  in  the  wake.  The 
interpolation  essentially  defines  the  distribution  of  vorticity  in  the  subdivided  interval. 
The  total  vorticity  in  the  subdivided  interval  still  remains  as  the  only  unknown  in  the 
wake.  Since  the  total  vorticity  was  unknown  in  the  original  2-D  formulation,  using 
Lagrange  interpolation  does  not  introduce  any  more  unknowns  to  the  boundary  value 
problem.  Instead,  the  interpolation  modifies  the  influence  functions  on  both  sides  of 
the  equation. 

The  sheet  vortex  strength,  multiplying  the  uniform  shed  wake  element  size  U8t, 
can  now  be  represented  by  Lagrange  polynomials  as  follows.  If  r  is  the  fractional 
time,  or  fractional  distance,  back  from  the  present  then 

A 

FaM    =     £  Lm(Tk)rusm  (A.S) 

m-0 

Lm(n)    =      n    Ih^rL  (A-9) 

where  n  represents  the  "product  of",  A  is  the  order  of  the  Lagrange  interpolation  and 
m  is  the  index  which  counts  the  un-subdivided  intervals  from  the  current  "zerotlT 
time  step. 

In  the  most  recent  time  step  element  (m  =  0),  the  discrete  subdivision  strengths 
can  be  obtained  by  applying  Lagrange  interpolation.  Care  must  be  maintained  to 
ensure  the  total  discrete  strengths  correctly  model  the  integrated  sheet  strength. 
Equation  A.S  and  A.  10  can  be  used  to  describe  the  unknown  subdivision  strengths  as 
a  combination  of  the  current  unknown  strength  and  the  known  strengths  of  the  prior 
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time  steps. 


Tsk  =  rA(rk) 


USt/Nj 
USt 


(A.10) 


For  example,  for  a  second  order  polynomial  the  Lagrange  interpolation  coefficients 
and  the  resulting  equation  describing  the  subdivision  strengths  are: 


Tsk  =  (PkTuSo  +  qkTuSl  +  rkrs2)US*/*f         with  k  =  1,  •  •  ■ ,  Nj 


Pk   = 


rl  -  3r,  +  2 


qk     =     -T2k+2rk 


rk     = 


n  -  n- 


(AM) 

(A.12) 
(A. 13) 
(A.14) 


It  should  be  clear  that  FuSi  =  Fsj^f  +  i  and  that  Tus2  =  Ts^f+2,  and  so  on.  which 
are  all  known  quantities  of  previous  solutions.  The  remaining  unknown  is  r"^0  which 
is  obtained  from  Kelvin's  theorem  written  as  equation  A. 4.  After  substitution  and 
collection  of  the  unknown  quantities  of  the  left  hand  side  of  the  equation,  the  formu- 
lation with  subdivision  can  be  written  as: 
N 


1     X/ 

-— y 

AT.   ^-^ 


*-tM^"» 


n-  -  tj 


■V 


(A.15) 


N 


/  fc= 


-r^-1  +  E[   II   — 


A  i  To  -  tj 


m=l    \j=Oj^m  Trn         TJ 


Nf+m 


A'u. 


k=Nf+l 


Or,  for  a  second  order  polynomial,  equation  A.15  simplifies  to  become: 


N 

E(r.) 

N 


B»'.i-  fEw-% 


N 


/  *=i 


=  -Vi 


(A.16) 


1        r  Nw 

""77"  y  [PkTn-i  +  qkTsNf+1  +  rkrsNf+2\  Wt,k  -     y    (rs)fcWi,* 

1   f  fc=l  k=Nf+l 

where  p,  q  and  r  are  defined  in  equations  A.12  through  A.14. 

Two-Dimensional  Results  with  a  Second  Order  Interpolation 

Like  figure  A-2,  figure  A-4  shows  a  arbitrary  time  step  for  two  cases.  The  first  case 
is  the  more  typical  single  vortex  representation  for  each  time  step  element.    The 
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Vorticity  Sheet  Strength  behavior  at  a 
random  timestep  in  a  sinusoidal  gust 


x/C  =  1.0isT.E. 


Figure  A-4:  Trailing  edge  behavior  of  the  vortex  sheet  strength. 

discontinuity  exists  at  the  trailing  edge  due  to  the  uneven  spacing  between  the  cosine 
spaced  foil  elements  and  the  constant  spaced  wake  elements.  The  second  case  shown 
in  figure  A-4  is  using  the  newly  developed  formulation  discussed  earlier.  The  new 
formulation  solves  the  boundary  value  problem  by  representing  the  first  wake  element 
more  like  a  sheet  as  opposed  to  a  single  concentrated  vortex.  Sheet  vorticity  further 
away  from  the  foil  a  sufficiently  represented  as  a  single  concentrated  vortex.  However, 
near  the  trailing  edge,  a  single  vortex  representation  of  the  time  step  element  causes 
erroneous  results  due  to  the  singular  nature  of  a  discrete  vortex  approaching  the 
trailing  edge  of  the  foil. 

The  major  advantage  to  this  formulation  is  that  it  allows  relatively  coarse  time 
steps  while  getting  the  accuracy  obtained  with  much  finer  time  steps.  Figures  A-5 
and  A-6  illustrate  this  advantage  by  examining  the  foil  vortex  sheet  strength  at  a 
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particular  x/C  location  as  the  gust  advances  over  the  foil.  In  these  figures,  the  first 
number  in  the  legend  is  the  time  step  element  size  relative  to  one  another.  The  sec- 
ond number  in  the  legend  is  the  number  of  elements  in  the  first  time  step  (A/).  For 
example,  figure  A-6  shows  a  case,  labeled  "2*,  1  element",  with  no  subdivisions  and  a 
very  small  time  step.  The  relatively  small  time  step  is  required  to  get  the  solution  to 
converge  to  the  correct  results  near  the  trailing  edge.  The  corresponding  subdivision 
case  to  comparison  against  is  the  "100*,  50  element".  With  subdivisions,  the  case 
converges  to  nearly  the  identical  results.  This  particular  subdivision  case  uses  50 
subdivisions  which  divide  the  "100*"  step  into  subdivisions  of  size  "2*".  Thus,  this 
case  closely  approximates  the  "2*,  1  element"  case.  The  advantage  is  that  the  extra 
50  subdivisions  are  solved  in  the  boundary  value  problem  and  are  not  required  to  be 
maintained  in  the  wake  past  the  first  time  step  element.  Additionally,  the  extra  sub- 
divisions do  not  increase  the  number  of  unknowns.  Instead,  the  subdivisions  simply 
modify  the  coefficients  of  the  simultaneous  equations  as  shown  in  equation  A. 15.  The 
new  formulation  dramatically  reduces  computational  effort  while  maintaining  similar 
accuracy. 

Conclusion 

The  new  formulation  dramatically  reduces  the  computation  effort  to  attain  the  same 
accuracy.  The  case  studied  herein  applies  a  second-order  interpolation;  higher  order 
interpolations  should  be  evaluated  to  determine  if  they  yield  even  more  computation 
gain  by  increasing  the  time  step  size  further.  Additionally,  this  method  should  be 
applied  to  a  three-dimensional  unsteady  code.  The  ultimate  use  for  this  method  is 
to  improve  real  life  calculations  in  a  three-dimensional  unsteady  method. 
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Wake  Spacing  Comparisions 

by  observing  the  foil  sheet  strength 

at  the  x/c  =  0.90  position 


c 

*— 
o> 
a> 


4  6  8 

Non-dimensional  time 


Figure  A-5:  Vortex  sheet  strength  behavior  at  x/C  =  0.9  on  the  foil  in  a  sinusoidal  gust.  Using  a 
non-dimensional  time  step  of  100*  with  Nj  =  10  vortices  in  the  first  interval  yields  similar  results 
as  a  time  step  of  10*  and  no  subdivision. 
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Wake  Spacing  Comparisions 
by  observing  the  foil  sheet  strength 
30  f~  at  the  x/C  =  0.98  position 


4  6  8 

Non-dimensional  time 


Figure  A-6:  Vortex  sheet  strength  behavior  at  x/C  —  0.98  on  the  foil  in  a  sinusoidal  gust.  The  new 
formulation  converges  to  similar  accuracy  while  using  significantly  larger  time  elements  in  the  wake. 
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Appendix  B 

Color  Figures 


PUF-14  Swept  Volume  Axial  Forces 
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Figure  B-l:  This  figure  shows  a  notional  inflow  with  the  corresponding  time-average  forces  in  the 
swept  volume  of  the  rotor  (see  section  6.3.3). 
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Underwater 
notional  body 
and  propeller 
at  30°  yaw  to 
port 

Color  contours 
represent  the 
computed 
axial  force 
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Figure  B-2:  Huang  body  1  coupled  lifting-surface/ RANS  solution  (see  section  8.4) 


L40 


J 

VT 

0.1519 

0.1060 

0.0600 

0.0141 

-0.0319 

-0.0778 

-0.1238 

-0.1697 

Figure  B-3:  This  figure  shows  the  volume  used  to  perform  the  stator  analysis.    This  case  is  at  4° 
drift  angle.  Tangential  velocity  is  shown  in  the  contour  (see  section  9.2) 
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Lateral  Force  (  F  ) 
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Figure  B-4:  Force  contours  on  the  body,  duct  and  both  blade-rows  (see  section  9.2) 
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